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l‘.l The theory of laminar instability is of great practical 

interest because it provides qualitative criteria for tlie impor- 
tant problem of laminar-turbtilent transition- It has been 
observed experimentally ^hatj laminar flow occurs at low 
Reynolds numbers and that in this range viscosity damps but 
any deviations from laminar flow. On the other hand, turbulent 
motion occurs at high Reynolds numbers. Laminar flow can be 
realised, if at all, only by excluding all possible disturbances 
One is led to suspect therefore that there might be more than 
one solution to the equations of fluid motion. The stability 
theory at least explains the noh occurrence of laminar flow 
undeh certain combinations of flow parameters. At present, the 
predictions of stability theory are mainly limited to the descri 
ption of the initial breakdown of the laminar flow due to 
infinitesimal disturbances. It is still a long way from 
explaining the regime of fully developed turbulent flow. 
Stability theory also predicts the qualitative effectiveness 
of the various flow parameters in promoting or suppressing 
stability. The fundamental aspects of the theory of laminar 
stability have been summarised in a monograph by Lin^. 
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ii2 ■ The classical theories of mechanics do not tale 

m . 

couple stresses into account. This essentially means 'that 
the mechanical interaction between portions of a body across 
a surface in it can be represented bv distributed forces 
only. In recent years much attention has been paid to polar 

O 

theories which take into account couple stresses. Stokes 
has considered one such polar theory for fluids and has 
studied the effects of couple stresses by solving some 
representative boundary value problems for fully developed 
viscometric flows. The results show that an interesting 
phenomenon involving a size dependent effect, which is not 
present in the classical nonpolar theory, comes in. 

In this thesis, an attempt has been made to study 
the effects of couple stresses on the stability of parallel 
flows. Numerical results have been obtained for plane 
Poiseuille flow. 

i,3 The mathematical problem of hydrodynamic stability 

can be formulated by using the given steady state solution 
of the equations of motion and superimposing on it a 
disturbance of suitable kind. This results in a set of 
non-linear. disturbance equations which govern the behaviour 
of the disturbance. If the disturbance ultimately decays to 
zero as time tends to infinity, the flow is said to be 
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stable. But if a disturbance results which is pe rmanently 
different from zero, the flow is said to be unstable. It 
does not follow that instability leads to turbulent motion, 
because another probably more complex form of the laminar 
motion may be the result. We assume that for small distur- 
bances the disturbance equations may be linearised. In 
general^ the resulting equations are homogeneous with 
homogeneous boundary conditions. This leads to a characte- 
ristic value problem with a parameter < 3 - . If all the 
characteristic values of cr have negative real parts, the 
motion is said to be stable with respect to infinitesimal 
disturbances. If some of the characteristic values of <7- 
have positive real parts, the flow is said to be unstable. 
Neutral stability is marked by zero values of the real part 
of (s— . The possible eigenvalues of cr depend, of course, on 
quantities such as flow speed, kinematic viscosity, thermal 
diffusivity and the wave length of the disturbance. 

A brief account of the development of polar 
theories is given in the next section. 

1.4 The classical theory of continua models properties 

of media in which a central force acts between particles. 

The inadequacy of this theory, for some problems, led to 
the consideration of theories of continuous media which take 
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irllb account rotational interaction amongst particles. This 
gave rise to the concept of couple stresses. It is well 
known that the action, at a point j of even a simple system 
of forces cannot be reduced to a single resultant force, but 
is in general equivalent to a force and a couple. Consequen- 
tly, to describe the mechanical interaction between a system 
of particles, it is necessary in general to consider forces 
and moments. 

The literature on couple stresses begins with the 

memoir of Cosserat and Gosserat . The Cosserats gave a 

systematic development of the mechanics of continuous media 

each of whose points have the six degrees of freedom of a 

rigid body. In classical elasticity, a material point has 

only three degrees Of freedom. 

A striking feature of Cosserats’ theory was the 

appearance of couple stresses in the equations of motion and 

the resulting asymmetry of the stress tensor. A modern 

derivation of the Cosserat equations was given by Truesdell 
A 5 

and Toupin^. Mindlin and Tiersten gave an extensive 
analysis of the infinitesimal motions of a Cosserat medium 
with constrained rotations. The mechanics. of elastic media 

with micro- structures has been discussed by Mindlin®. 

. 7 

Eringen and his coworkers developed a theory for 
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microscopic materials whicii exhibit Certain microscopic 
effects- arising from the local stihictiire and micromotions of 

. g ,r 

the media; Eringen also simplified this general theory to 
the micro polar theory and solved some boundary value 
problems. Upadhyaya ’ studied the effects of couple 
stresses in linear viscoelasticity; 

1,5 Many theories of polar fluids have been developed, 

A polar fluid is a fluid in which couple stresses are present. 
In the general polar fluid theories, the angular velocity 
of a fluid particle is allowed to differ from the average 
angular velocity of its neighbourhood (which is called the 
vorticity) by introducing an additional independent angular 
velocity (or internal spin) vector field. By constraining 
the rotation of a fluid particle to equal the local rotation 
of the medium, the couple stress theory results as a special 
case of the general polar fluid theory. 

The theory of polar fluids is analogous to the theory 
of polar elasticity in the same way as the classical theory 
of elasticity is analogous to the classical theory of fluids. 
Grad^^(l952) first introduced the constitutive equations for 
a polar fluid. These constitutive equations relate the skew 
symmetric part of the usual stress tensor to the relative 
angular velocity of a fluid particle (with respect to its 
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neighbourhood), and the couple stress tensor to the gradient 

of the total or particle angular velocity. The Newtonian 

law of viscosity is preserved by assuming that the usual 

relationship between the s 3 nmnetric part of the stress tensor 

and the rate of deformation tensor holds. Aero, Bulygin and 
1 P 

Kuvshinkii introduced the velocity and total angular 

velocity fields and postulated a dissipation function that 

l3 

led to Grad’s constitutive equations. Cowin (1962) used 

a Cosserat continuum in which a rigid triad of vectors is 

associated with each particle. The total or particle 

angular velocity is then represented by the angular velocity 

of the rigid triad and leads -to the same constitutive equa- 
14 

tions. Eringen (1964) developed the theory for micro- 
fluids. By specialising the micro-fluid theory to the 

15 

theory of micro-polar fluids, Eringen obtained constitutive 
equations similar to those of Grad. Dahler developed the 
same constitutive equations from a statistical mechanics 
approach. 

Stokes obtained a linear theory for fluids with 
couple stresses parallel to the couple stress theory in 
linear elasticity proposed by Mindlin and Tiersten, The 
symmetric part of stress tensor and the rate of deformation 
tensor follow the same relationship as in the non-polar 
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case. The couple stress tensor is assumed to be proportional 

to the curvature twist tensor which is defined to be vorticity 

gradient. This theory is a special case of the general polar 

theory, obtained by constraining the rotation of a point in 

a Cosserat continuum to be equal to the local rotation of 

the medium. The antisymmetric part of the stress tensor and 

the trace of the couple stress tensor are left undetermined 

by these constitutive equations but can be determined from 

17 

the boundary conditions. Condiff and Dahler and 
X3 

Pennigton"^ have considered plane Couette flow, plane 

Poiseuille flow, Couette flow between two concentric cylinders 

l3 

and Poiseuille flow in a pipe. Cowin considered plane 

Poiseuille flow, Eringen and many others have considered the 

flow of micro-polar fluids in different flow situations, 

Stokes has solved Couette and Poiseuille flows between 

parallel plates, flow between concentric cylinders with and 

without a toroidal pressure’ gradient and Poiseuille flow in 

a pipe. Stokes ’ also considered the effects of couple 

stresses in fluids on hydromagnetic flow in a channel and on 

2l 

heat transfer. Dep has studied the equations for fluid 

boundary layers with couple stresses. Dep's analysis is 

based on the equations given by Aero, Bulygin and Kuvshirikii. 
22 

Usmani considered the time dependent effects in fluids 



8 


with couple stresses, using the theory proposed by Stokes. 

In this thesis an attempt has been made to study 
the effects of couple stresses on the stability of parallel 
flows. It has been assumed that the fluid is incompressible 
and that the body forces and body moments are absent. The 
equations of motion are fourth order partial differential 
equations. The method of normal modes has been used to 
study the variation of disturbances with time. 

The disturbance equation comes out to be a sixth 
order ordinary linear differential equation subject to homo- 
geneous boundary conditions. This leads to an eigenvalue 
problem. The stability of parallel flows has been analysed 
for three-dimensional and two-dimensional disturbances. It 
has been observed that the minimum critical Reynolds number 
is obtained in the case of two-dimensional disturbances. 
Keeping this in view, stability of plane Poiseuille flow 
has been studied numerically. The initial value technique 
suggested by Nachtshiem^^as been used to solve the charac- 
teristic value problem. The integration is done by a step- 
^^hegration method. To have better control on 
truncation error, the fifth order Milne predictor-corrector 
integration technique is used. Double precision arithmetic 
(16 significant digits) has been used to control round-off 
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.errors. The effects of couple stresses on the Reynolds 
number, on the disturbance wave number and on the stability 
factor have been studied. Neutral stability curves for 
different values of couple stress parameter have been 
obtained. The minimum critical Reynolds number in each 
case has been obtained. 



2 . HYDRODYNAMIC STABILITY WITH GOUPLE STKS SSBS 

p 

Following Stokes , the theory of couple stresses 
in fluids is briefly discussed in this section. 

2.1 KINEMATICS OF FLOW : 

Let UjL be the components of the velocity field. 
Then, the rate of deformation tensor Dj_:^ and the vorticity 
tensor are given by, 


^13 

_ 1 

“ 2 



(2.1.1) 

w. . 

13 

1 

2 

- 


(2.1.2) 


respectively. 

The vorticity vector is defined as half the curl 
of the velocity vector, so that 



^i 

1 

= 2 ®iis ^s^r 

(2.1.3) 

The 

vorticity 

vector and the vorticity tensor are 

related by 

the 

equations 




%3 

eijr 

(2.1.4) 



1 

(2.1.5) 

and 

Oi 

= — e. W 

2 irs ^rs 


The curvature-twist rate tensor is defined to be the 
gradient of the vorticity field; thus, 



( 


11 


The diagonal components of are a measure oi the rate of 
twist of the material per unit length about the three axes; 
On the other hand, the off diagonal elements of are a 
measure of the curvatures induced per unit time in the 
various planes. 


liquations (2.1^3) and (211.6) give, 

= 0 _ (2.1.7) 

That is, the average twist rate per unit length of three 
mutually perpendicular line elements is 2ero. 


2.2 CONSTITUTIVE EQUATIONS: 

In the non-polar case there is a close similarity 

between the elastic and fluid constitutive equations. 

Following Mindlin and Tiersten' s linear constitutive equa- 

2 

tions for polar elasticity Stokes proposed the following 
constitutive equations for fluids, 


T®. = - P<?i3 Drr ^ (2.2.1) 


+ 4T|'E^j_ (2,2.2) 

where the stress tensor Tj_^ and the couple stress tensor 
are given by the equations, 


n m , . 


^i = 


(2.2.3) 


where n. is the unit normal to the surface on which t^ and 


mj_ act. 
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It is clear from equations (2,2.1) and (2.2.2) that, 
for couple stresses to exist, the vofticity gradient field 
must he non zero. The dimensions of the material constants 
A and are those of viscosity (namely M/LT) , whereas the 
dimensions .of '^1 and are those of momentum (namely ML/T) . 
The ratio has dimensions of length square. This 

material constant is denoted by where 


i= ' (2.2.4) 


For incompressible fluids u„ _ = 0, and in the absence of 

r,r 

body forces and body moments, Cauchy's law of motion 
governing the balance of linear momentum is given by. 





rr 


n 


u. 


i,rrss 


2.3 BOUNDARY CONDITIONS'. 


(2.2.5) 


In dealing with this polar theory six boundary 
condit‘’’ons are necessary. It seems quite reasonable to 
assume the no- slip condition at the boundary as in the non- 
polar case. The formulation of additional boundary condi- 
tions is restricted to certain idealised cases only, 
because the mechanism of interaction at the boundary is 
not understood clearly. In the general theory of polar 
fluids the most commonly used boundary conditions are for 
limiting case wherein a solid surface interacts so strongly 
with the neighbouring fluid, that a particle of fluid does 
not turn over relative to the surface, and therefore, its 
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angular velocity is equal to th.e angular velocity of the 
surface. This makes use of the property of adherence of a 
fluid to a solid surface. In the constrained polar theory, 
this condition becomes trivial, as the particle and the 
local angular velocities are assumed to he same throughout 
the flow field. 

To prescribe additional independent boundary 
conditions in this case, use is made of the other idealised 
situation wherein the solid surface allows rotational slip 
of a fluid particle in such a way that the couple stress 
vector is zero at the surface. This leads to boundary 
condition ; 

^i ” ^^3 ^ji j surface ~ 

where n^ is the unit normal to the surface. 

p 

Following Stokes the condition of no couple 
stresses at the solid boundaries is used here._ 


2,4 GOVERNING EQUATIONS IN NON-DIMENSIONAL FORM ; 

For incompressible fluids, in the absence of body 
forces and body moments, the velocity components u^ and the 
pressure p* satisfy the equation of continuity, 

u^ 0 (2,4.1) 

i , 1 

and the equations of motion, given by equations .(2 .2.5) , 
namely, 



(2.4.2) 
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The asterisk indicates that the quantity concerned is 
dimensional. The superposed' dot indicates the derivatives 
with respect to time following a particle, and ^ is the 
mass density of the material. 

It is convenient to work entirely in terms of 
non-dimensional quantities. For this purpose, following 
characteristic quantities are introduced ; 

L = characteristic length (length scale) 

'V* = characteristic velocity (velocity scale) 

L*/v* = characteristic time (time scale) 

_ characteristic pressure (pressure scale) 
The non-dimensional quantities may be obtained 
from dimensional quantities by making use of the simple 
relation % 


, ... Dimensional quantity 

Non-dimensional quantity = s^ale — 

The non-dimensional variables are defj.ned by, 
t = 

Xi = 

u^ = Uj_*/V* 

p = p*/ 

The equations (2.3.1) and (2.3.2) may be written 
in non-dimensional form as, 


u. 


1,1 


0 


(2.4.3) 
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and 


^ ^i,r = -P,i + I ^i,rr “ ^ ^i,rrss 


(2.4.4) 


The non-dimensional parameters R and K are given by. 


R = V* L* ^ 


(2.4.5) 


K 


L * 2/^2 


Here R resembles the Reynolds rnimber in non-polar theory 
for fluids. 

At solid bon.ndaries with prescribed velocities Sj_ 
and couple stress tensor Cj_^ the conditions 

Ui = Si 


and . (2.4.6) 

"1 = ’'3 °i± 

must be satisfied. 

If the boundary velocities and couple stress 
Gj_^ are independent of time, we may have a steady motion 
given by, 


% = Hi 

P = P Cxj^) 

njMji = nj Mji(xj,) 


(2.4.7) 


Satisfying equations (2.3,3) and (2.3.4) in' the form, 

1 

R “i,rr bk “i,rrss 


U4 u. . 


- p . + i u_. ~ Uj 


''\,k 


= 0 


(2.4.8) 

(2.4.9) 
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together with the boundary conditions, 

Ui = S;J_ 

(2.4.10) 

5. M AT HEMIT IG AL K) EMULATION OF THE STABILITY PROBLEM FOR 
AN INCOMPRESSIBLE FLUID ; 

To study the effects of small disturbances on the 

steady motion, we seek solutions of equations (2.4.3) and 

(2.4.4) in the form 

^i “ ^i ^ * 

(2.5.1) 

P = P + g; ^ + 

where £ is a constant parameter and u^^\ ..... 

( o\ 

p^ , ..... are functions of position and time. We demand 
that equations (2.5.1) satisfy equations (2.4.3), (2,4,4) 
and (2.4.6) for all values of £ in the range 
Formal substitution from equation (2.5.1) into equations 
(2.4.3), (2.4.4) and (2.4.6) gives a sequence of sets of 
four equations, each set corresponding to a definite powder 
of £ . The set corresponding to £^ is given by equations 
(2,4.8) and (2.4.9) along with the boundary conditions 
given by equation (2.4.10). The set corresponding to £,-^ is 

- ti) (1) - (1) 1 (1) 

= + U. U. . + U. U . . = - p .. + p U- 

5t 3 1,3 3 1,3 ^ 


< 2 - 6 - 2 ) 
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and the continuity equation is ' 

= 0 (2.5.3) 

The boundary conditions for this set are 


= 0 


n. = 0 


on solid boundaries. 


A complete treatment would require consideration 
not only of equations (2.5.2) but also of the equations 
corresponding to higher powers of £ , together with the 
consideration of the convergence of equations (2.5.1) and 
a 3'iistif ication of the term-by-term differentiation. The 
problem of hydrodynamic stability reduces to the problem of 
solving a system of non-linear (quasi-linear) partial 
differential equations. 

In the usual approaches to the problem, the 
mathematical formulation is cast in a different way. ¥e 
assume that, for small disturbances the equations may be 
linearised! that is, we neglect quadratic or higher order 
terms in the disturbances and their derivatives. Thus, 
the linearised theory corresponds to equation (2.5.2). The 
linearised theory of stabilitv of laminar flows decomposes 
the motion into a mean flow (whose stability constitutes 
the sub.iect of the investigation) and into a disturbance 
superimposed on it. 
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Let the mean flow be given by and pressure by p. 
The corresponding quantities for the non-steady disturbance 
will be denoted by u^ and p' respectively. The resulting 
motion is given by, 

Uj_ = Uj_ + u^ (2.5.5) 

P = P P ^ (2,5 .6) 

The couple stress tensor may be written as, 

^ 2 . 6 . 7 ) 

It is assumed that the disturbance quantities are small 
compared with. the corresponding quantities of the mean flow. 

The investigation of the stability of such a 
disturbed flow can be carried out with the aid of either of 
two different methods. The first method (Energy method) 
consists merely of calculating the variation of the energy 
of the disturbance with time. Conclusions are then based on 
whether the energy decreases (the flow is stable) or 
increases (the flow is then unstable) with time. The theory 
admits an arbitrary form of the superimposed motion and 
demands only that it should be compatible with the equation 
of continuity. In general, the energy method has not 
proved successf'^l* The second method accepts only flows 
in which the steady state and total velocity fields satisfy 
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the equations of motion. This is the method of saall 
disturbances and has proved successful i 

It is assumed that the mean flow given by equation 
(2.4.7) is a solution of equation (2'.4v4) and that the 
resultant motion, governed by equations (2,5.5), (2^5,6) and 
(i2|5,7), also satisfy this equation. Substituting for the 
resultant motion in the equation (2.4.4) and linearising, we 
get 


dui . - ..t 


i 4,r + 4 "l,r = ' 


a* 


The continuity equation becomes 


u. . = 0 

1,1 


i-ul 

RK i,rrss 


(2.5.8) 


(2.5.9) 


The boundary conditions reduce to. 


ui = 0 


n^ M’i = 0 


(2.5.10) 


A study of the linear system of equations (2.5.8) shows that 
the time appears only through the derivatives with respect to 
time t and hence solutions containing an exponential time 
factor e'^ may be expected. Since the original boundary 


conditions are already satisfied by the steady motion, the 
disturbance solutions must vanish at the boundaries. This 



20 


■will be possible only for special values of cr . Thus we have 
a characteristic value problem with 5- as a parameter, u^j 
p* and may be referred to as the modes of oscillations. 
Usually, by slightly changing some of the flow parameters, 
such as the Reynolds number, one finds slightly different 
eigen solutions for u^, p^ and by allowing a— to be 
complex. If the complex eigenvalue s-has a positive real 
part, the solutions u^, p' and are amplified because 
they tend to infinity as t increases and instability is 
indicated. If r' has a negative real part, the solution is 
’damped* and the steady flow is stable. The transition 
from stability to instability, in such cases, is marked by 
the vanishing of the real part of o- . It may -be noted, 
however, that a degenerate case of the transition from stabi- 
lity to instability occurs when the real and imaginary parts 
of <j~ Vanish simultaneously. The neutral mode then corres- 
ponds to a ‘secondary flow' rather than to a 'steady 
oscillation’. In certain problems, it is assumed that the 
imaginary part of tr- always vanishes if the motion becomes 
unstable. This can be shown explicitly for some problems. 

The secondary flow simply diverges with time after the 
stability boundary is crossed. In this case, one locates 
the stability boundary by finding the conditions for the 
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existence of steady secondary flow. 

We shall formulate the stability problem for 
parallel flows in next section. 



3* STA BI LITY OF PARALLE L FLOWg 

» 

3.1 fomulation of stability problem for parallel FLOVJS 

Let the fluid occupY the fegibn hatv/een two infinite 
parallel plate which are either stationary or moving 
parallel to each other with uniform speeds. 

Take a cartesian coordinate system with the x axis 
parallel to the direction of motion and the y and z direc- 
tions perpendicular to it. The plates are situated at 
y = y^ and y = 

The mean velocity field is chosen as? 

U = u (y), 0, 0 (3,1.1) 

The disturhance field is defined as, 

U’ = (u* , v' , w’) (3.1.2) 

p— p+p* (3.1,3) 

where p = P (x) 

The disturbance quantities are functions of x, y, z and time t. 
The governing equations (2.5.8) and (2,5,9), for the flow 
situation governed by equations (3,1.1) and (3.1.2), reduce 


to 
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4- -y* 

dt ax ay 


^ + 
oX 






y + 
at 


av* 

ax 


dp 

ay 


+ 


e'^ 


1 


y - 


4 , 

V ^ 


aw ' , — dw' 

at ^ ^ ax 



1 

R 



w 


1 

RK 



w 


(3.1.4) 


where 

and 

The boundary conditions for the system of equations (3.1.4) 
become 


u’ = 0 ■ 

v' = 0 

w' = 0 at y = y^ and y^ (3.1.5) 

n^-M 1 . = 0 
3 11 

The precise discussion of this problem presents certain ' 
difficulties. Physically speaking, no experiments are 
performed x\rith an apparatus of infinite dimensions. Mathe- 
matically, it is possible to circumvent the difficulty by 
limiting the discussion to disturbances which are spatially 
periodic in the directions in which the fluid extends to 
infinity. This is suggested by the nature of equations 
(3.1.4). Since all the three variables x, z, t are cyclic 
and since the coefficients of equations (3.1.4) depend on y 
only, the system admits solutions which are exponential 
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functions in x and z, as well as in.t. If the solution is 
to he hounded for x and z both tending to + oo and -ooj the 
corresponding exponents must he pure imaginary. Thus, each 
component of the disturbance is the real part of an expression 
of the form 


q/ = q'(y) exp |i (o(x + jBz) - i.x'ct j. 
- lore = (T? L*/V* 


(3.1.6) 


where c< and yS are real. This type of representation is 
generally referred to as the ’method of normal modes’. The 
three dimensional disturbances given by equations (3.1»6), 
having all the three components may be interpreted as 
oblique plane waves travelling in a direction making an 
angle uf tan”^( ) with the main flqV direction. 


Since the directions are arbitrary, it is convenient 
to choose them in such a way that, the x axis coincides with 
the direction of the wave front. 

The complex non-dimensional characteristic value a~ 
has the value (T = - icKc, in equation (3.1.6). The real 
Wave numbers 'X and /3 are given by, 

<X = 2 

/? = 2 7r/A 2 

Where and are non-dimensional wave lengths in the x 
and z directions respectively. 
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3.2 STABILITY CRITEIIIA 

In the expression given hy equation (3.1,6) the 
complex quantity c may he written as 

= = Oj, + 1 Oi 

where the real part represents the wave velocity of the 
disturbance and the imaginary part c^, the amplification or 
the damping of oscillations with time. The eigenvalue g— 
may be written as 

d~ = CT"^ + i cr 

= - i <KCy 

As discussed in section 2,5, the stability criteria can be 
based on the real part of as, 

(1) 0^ > 0 the flow is unstable, 

(2) ci = 0 neutral stability, 

and (3) 4 0 the flow is stable. 

3.3 TWO DIMENSIONAL AND THREE DIMENSIONAL DISTURBANCES 

From equations (3.1.6), the velocity and pressure 
terms may be written as 

|u’, v', w’ , p’j = Real part ) u(y), v(y), w(y), 

p(y) j I exp L i( x+ /3 z) ~i !< ctj I 

(3.3.1) 
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The equations (3.1.4) qan be simplified by substituting the 
values from equations (3.3.1). The simplified equations come 
out to be, “ - 




— 




• L 

Li‘ 

(u - 

c) 

A 

U = 

1 

. — . .\ A 

(Du) V + 1 ^ p 


- - io< 
R 

(I - 

c) 

A 

V = 

A 

Dp 

(3.3.2) 

l-io< 

__R 

(u - 

c) 

A 

w = 

i^p 



. ■ A A A . ^ 

i( c<u + pv) + Dv = 0 
where the operators are defined by, 

D - d/dy 

Is - (oC^ 4-^^) (3.3.3) 

and L = L (1 - L/K) 

The above equations hold for three dimensional 
disturbances. Following Squire (1633), we shall now show 
that the problem of the- three dimensional disturbances is 
actually equivalent to a two dimensional problem at a lower 
Reynolds number. In fact, if we introduce the following 
transformations, 

squ = c<u + (3 w 

V = V 
A 

p R = pR 
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g?(R “ c<R 
c = c 

5-2 =o<2+/j2 

The equations (3.3.2) may be written as 


(3.3.4) 


L - io< R (u - c) 


ic^ R (u - c) 


^ - R (Du) V + i .^RP 


V = R D p 


where 


ic<u + D V =0 
^ (D^ -^^) (1 - (D^ -22)/k) 


(3.3.5) 


This system of equations has the same structure as eqiua- 
tions (3.3.2) with w = 0,/o= 0. That is, equations (3.3.5) 
correspond to a two dimensional disturbance. The boundary 
conditions obviously reduce to, 

= V = 0 

^ at y = y and y (3.3.6) 

n. M.. = 0 1 2 

3 31 

Thus, each three dimensional problem is equivalent to a two 
dimensional problem. Hence it is sufficient to solve 
two dimensional problems only. Three dimensional solutions 
may then be obtained by using the transformations indicated 
by equations (3.3.4). In fact equations (3.3.4) show that 
the equivalent two dimensional problem is associated with 
a lower -Reynolds number, since U >oC. Thus, the minimum 
critical Reynolds number is given directly by a two 
dimensional analysis. 
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Another method of arriving at the above result is 
physically more suggestive. The three dimensional .distur- 
bance is essentially an oblique wave, how, if the x axis is 
chosen to be such that /3 = 0, we see that the equation for 
w(y) is uncoupled from the rest of the equations and u(y) 
and v(y) can be solved independently of w(y). As the 
direction of the obln.que wave is changed, the shape of the 
velocity profile u(y) remains unchanged, however, the 
magnitude is given by the pro;iection of the actual velocity 
vector on that direction. Thus the wave propogation in the 
direction of the velocity vector U(y) corresponds to the 
highest effectj.ve velocity or Reynolds number. 

Usually a laminar flow is always stable at 
sufficiently low Revnolds numbers because of viscous dissi- 
pation. A finite critical Reynolds number may be reached 
when at least some self- sustained oscillations become possible. 
In the present case, it then follows that, the critical 
Reynolds number for an oblique wave must be higher thc-n 
that for the wave propagating in the direction of the 
steady mean flow. Because of the uncoupling, one needs to 
consider two dj.mensj.onal disturbances only. This is a very 
Important result, as it Justifies the consideration of two- 
dimensional disturbances only. The uncoupling of the 
lateral oscillations is possible only when the curvature of 
the surface is negligible. The stability boundary for the 
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oblique wave is obtained by shifting the stability boundary 
for the two-dimensional disturbances towards higher Reynolds 
numbers at each wave number o<_ i As a result, however, 
certain combinations of o< and^- R may become unstable for the 
oblique wave although stable for the two dimensional 
disturbances; 

3.4 STABILITY OF PARALLEL FLOWS WITH RESPECT TO 

longitudinal-vortex disturbances and transverse- 
vortex DISTURBANCES; 

Disregarding the lateral disturbances w', the 
disturbances u’ , v^ may be visualised as being equivalent 
to a disturbance vorticity distribution which is periodic 
in X and is in the direction of the z axis only. We may 
refer to this type of disturbances as 'Transverse-vortex 
disturbances'. 

On the other hand, suppose we consider the alter- 
native extreme of setting o(= 0. Then the roles of u' and w’ 
are exchanged (while dealing in cartesian coordinates only) 
and nothing new follows. However, the disturbances v' and w' 
are solved independently of u’ . These disturbances may be 
defined as 'Longitudinal-vortex disturbances' , the vorticity 
component being parallel to the x axis. Although the longi- 
tudinal-vortex disturbance is of trivial importance for 
parallel flows, over a flat surface with steady state 
velocity in the x direction only, it turns out to be of 
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great importance in case there is a centrifugal force due to 
the curvature of the surface in the x direction or in case 
other body forces such those due to thermal buoyancy exist. 


3.5 IIIALYSIS OF TWO DIMENSIONAL DISTURBANCES: 

For a parallel flow with velocity u (y) under 
t ran averse -vortex disturbances, it is only necessary to 
consider two dimensional disturbances u' and v’ . The 
equations of motion are obtained from equations (3.3.2) and 
(3.3.3) by putting |3 = 0 and w' = 0. The operators given 
by equation (3.3.3) reduce to, 

D = d/dy 

L = (d2-o(^) (3.5.1) 

L = L(l-L/K) 

and the equations (3.3.2) reduce to, 


. R 

L 

R 


i cs< (u-c) 


icX. (u - c) 


/\ — /V , 

U = (Du)v 4- 1D(P 


A 

Y 


/V 

Dp 


By eliminating the pressure term and by making use of the 
continuity equation, the above system of equations reduces 
to , 


i i Lv ! = i(KR 


R 

Subject to boundary conditions: 


■_ , 2_“U 

(u-c) L - D u V 


(3.5.2) 


u = V = mJi = 0 at y = y^ and 
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In "the limiting case of vanishing couple .stresses the 
non-dimensional parameter K— cxc • so that equation (3.5.2) 
takes. the form, 


(d 2 V = io< R I (u-c)(D^. - D<,^) V - (D^)v 1 

f. I 

(3.5.3) 


Subject to the boundary conditions 

n = v = 0 aty = y^ and 

Equation (3.5.3) is the ’Orr-Sommerfeld' equation. 
Equation (3.5,2) is essentially an equation for 
vorticity. A general tvro dimensional motion can be speci- 
fied by a stream function 

(x, y, t) such that, 


u = 

V = - ^x 

The vorticity is then given by, 


(3.5.4) 


r = ^ ^ V 

^ b7 " dx 

= A-T 

where A = 

In terms of the vorticity , equations (2.5.8) reduce to, 


^ C ^ 5 T ht 
at ay ax 


LX 

d X 


ay 




(3.5.6) 
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To study the effects of small disturbances on the 
steady flow given by Y (3c,y) we make the substitutions 


y (x,y,t) =■^y(x,y) +V*(x,y,t) 

and (x,y,t) =^(x,y) + ^’(x,y,t) 

5-n equation (3.5.5). Linearising the resulting 


(3.5.6) 


equations, 


we get, 


iii 

at 


+ u 



sx ^y 


i A ^ 

R ^ 


i-zk At; 

RK 


(3.5.7) 


For parallel flows y''(x,y) is a cubic in y and is indepen- 
dent of X. 

The different j. adequation for 'S,’ is linear and has 
coefficients independent of x and t. Consequently, we may 
expect a solution of the form, 

V (x,y,t) = 0 (y) (3.6.8) 

The disturbance vorticity is then obtained as, 


C = (d2-o< 2) 0 (y) 

and equation (3.5.7) takes the form 

(d2_o<;2)2 (i_(d 2_ 0 _ jju-c) (D^. ^ 

- (D^ u) (3.5.9) 

The appropriate boundary conditions are 

u’ =0 = 0 

v' = D0 = 0 

m' = (D^ - C<^)D0 = 0 
yz 


at y = y^ and y^ 
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The aniplihude S of the vorticlty function is given hy, 

S = (D^ - 0^2) 0 (3.5.10) 

Therefore, the boundary conditions can be -written as 

0 = 0 

W = 0 at y = y^ and y^ (3.5.11) 

DS = 0 

The equation of motion (3.5.9) subject to bo-undary conditions 
given by equations (3.5.11) results in a characteristic value 
problem of the type 

F (CK, R, K, c) = 0 (3.5.12) 

For each combination ofo(y R and K there is a characteristic 
value c which is in general complex. For fixed values of K, 
the condition c^ = 0 leads to a relation between cKand R, 
that is, to a curve in theo<^, R plane. This curve is usually 
referred to as the neutral stability curve. 

The stability of plane Poiseuille flow will be 
discussed in the next section. 



4. STABILITY OF PLAjSfE POIBEUILLE FLOW 


4.1 REVIM OF LITERATURE ON STABILITY OF PLANE POISEUILLE 
FLOW ; 

The class of strictly parallel flows is limited, 
since u(y) can at most he quadratic in y. This includes, 
however, two important special cases ; 

(1) Plane Couette flow ; 

n(y) = y . (-i ^ y ^ l) 

where V* = The velocity of the upper plate 
L = Half the channel width 

(2) Plane Poiseuille flow •, . 

u(y) = 1 - y^ (-l^y 4l) 

where V* = The maximum velocity at the centre of 
the channel 

L* = Half the channel width 

In plane Couette flow, the fluid motion is a simple 
shear caused hy the relative motion of the two plates. This 
is obviously the simplest laminar motion conceivable, and one 
night expect a simple stability problem. However, no conclu- 
sive answer has yet been arrived at for this problem. All 
existing investigations tend to show that the flow is stable, 
'n the case of plane Poiseuille flow an exact solution of the 
)roblem, even in a formal sense, is not possible. This flow 
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therefore provides one of the "basic problems to be investigated 
for stability. 

Historically the study of laminar flows, such as 
classical plane Couette flow and plain Poiseuille flow was 
started with the consideration of a transverse-vortex distur- 
bance by Rayleigh^^ followed notably by Heisenberg^"^, 
Tollimien^^ and Lin^®. 

The stability of plane Poiseuille flow with respect 
to two-dimensional disturbances has been studied by many 
authors. A considerable amount, of controversy exists because 
of the contradictory conclusions which have been arrived at. 
Helsenberg^"^ concluded that plane Poiseuille flow becomes 
unstable at a sufficiently high Reynolds number, but did not 
arrive at a critical value beyond which instability begins. 
Lin^ obtained the minimum critical Reynolds number of 5300, 
based on the maximum velocity at the centre of the channel 
and its half width. Calculations of curves of constant 
were made by Shen^"^ by a perturbation from Lin’s neutral curve. 
Both Heisenberg and Lin used. as 3 nnptotic series. Mel£S3m28(i946) 
used a similar method and obtained a minimum critical Reynolds 
number of 6800 after certain approximations. A different 
method was used by Pekeris^^, who concluded that plane 
Poiseuille flow is stable at all Re imolds numbers. It is 
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believed that this conclusion arises from the fact that 

Pekeris’s series is essentially valid only for stable cases. 

The disagreement of Pekeris's result with Lin's calculations 

led, Von Neumann to suggest a direct numerical calculation 
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which Was made by Thomas . In 1953, Thomas, using the Finite- 
difference technique, arrived at the minimum critical Reynolds 
number of 5780. The Orr-Sommerfeld equation of hydrodynamic 
stability theory has presented a challenge to both asymptotic 
and numerical analyses for almost half a century. In spite of 
its appearance as a simple linear eigenvalue problem, its 
structure has required development of special asymptotic 
methods as well as numerical techniques. The asymptotic 
methods are described in detail by Lin^. The asymptotic 
theory of the Orr-Sommerfeld equation has been greatly advanced 
in recent years. However, it is incomplete in many respects. 

In much of the recent work on this aspect of the problem (Lin 
and Rabenstein^l(l960) , ' Sibuya^^(l963)), the major emphasis 
has been on the uniformity of the approximations achieved. 

It has been shown that uniform approximations to the solution 
of the Orr-Sommerfeld equation can be obtained in terms of 
the solutions of a simpler comparison equation. This equation, 
however, is also of the fourth order. All of this work has 
been restricted to flows with one critical point for which 
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the comparison equation can te solved by the method of 
♦ 

Laplace integrals. Even for such flows, the theory is quite 

complicated aind its extension, for example, to flows vrith 

more than one critical point, though perhaps not formally 

difficult, would clearly lead to further complications. ' In 
33 

1966, Graehel applied the technique of matched expansions 
to the classical inner and outer solutions to obtain eigen- 
value relations. The technique of Multiple- scale was used by 
34 

Kuentam to obtain the asymptotic solution of the Orr- 
•Sommerfeld equation. The stability of steady and time depen- 
dent plane Poj.seuille flow has been studied by Grosch^^ by 
making use of the method of expansions in a set of orthogonal 
functions. Much work has also been done in obtaining numerical 
solutions of the stability problem. The direct numerical 
calculations made by Thomas were quite lengthy and, according 
to Lin^, the limited amount of work required two weeks of 
machine time on a high speed e.lectronic calculator. Nachtsheim^^ 
successfully applied the initial value technique in solving the 
stability of plane Poiseuille flow. Nachtsheim obtained the 
minimum critical Reynolds number of 5769. The agreement of the 
results of Nachtsheim with those of Thomas is very good. 

Kaplan^'^(l964) also applied the Finite-difference technique for 

38 

calculating the eigenvalues. Lee and Reynolds suggested two 
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numerical techniques to solve the 0 yr- Somme rf eld equation. 

The first is a variational approach which allows the eigen- 
value problem to be reduced to an algebric problem of matrix 
eigenvalue determination. The second scheme involves numeri- 
cal integration of the differential equations. This scheme 
is based on the method suggested by Kaplan with some modifi- 
cations. Potter , making use of ?ce 3 molds technique, 
studied the finite amplitude stability of parallel flows. 

In the following sections, the numerical scheme for 
solving the stability of plane Poiseuille flow with couple 
stresses has been developed. Bie method suggested by 
Nachtsheim is the basis of the present work, 

4.2 FOEMULATION OF THE PROBLEM ; 

Consider the flow due to a pressure gradient between 
two parallel plates at a distance 2h. Let the x^ axis 
coincides with the lower plate. The axis X 2 is normal to x^ 
axis and the Xg axis represents the lateral direction. The 
governing equation of motion in a dimensional form is, 

= - P,i * 

Body forces and body moments are assumed to be absent. For 
steady one dimensional flow of sm incompressible fluid between 
two parallel plates, the velocity field is given by 
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Ui = u^Cxg) 

Ug = 0 

U 3 = 0 

Equation (4.2.1) then reduces to, 

1 d^u. . dp 

ax2^ ""~L ^ (4.2.2) 

Where p = p(x^) 

p 2 _ 5 ■^being a material constant. 

The general solution of equation (4.2.2) is 

h = ^o * H ^2 +2^ ^2 Bi=°sh(Xg/£ ) 

+ Eg sinh (X 2 /^ ) (4.2.3) 


The. only non-zero component of the vorticity, defined by 


6 :>.? = 1/2 eirs u^ _ , is given by 

-L S ^ X 


to O - “ 

2 2 


^^1 

52:2 


the couple stress tensor is given by 

“13 = 


where K. . =63. . 

•^o tJ ^ 1 

the only non-zero component of interest is? 
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M = 4 7i K" 

23 ^ 1^23 


2 r' i 


3- (^) + !i 

H ^ 




B2 




cosh_ (Xg/ ^ ) 


- (— ^) slim (x 2/£ ) 
besides this, M^g is also non-zero. 

The appropriate boundary conditions are, 


= 0 


M 


23 


0 


at Xg - 0 and Xg = 2h 


(4.2.4) 


(4.2.5) 


Solving equations (4.2,3) and (4,2,4) with boundary conditions 
given by equation (4.2,5), we obtain. 


Ul(Xg) = (~^)| 

12 2 - (^x^ [ 


2 2 

“ 2 (xg/h) + (Xg/h) + ~ h - cosh 

3. L 


X2 '■) 

(xg/ ^ ) + tanh a. sinh (— -) { 


(4.2,6) 


The maximum velocity occurs at the centre of the channel, 
Xg = h and is given by, 


"'^Imax 



(!£-) 

3^1 



^ (l - l/cosh a) 

3 


where a = h/o (4.2.7) 

Equation (4.2,6) is non-dimen sionalised by taking li^jnax; 
the characteristic velocity and half the channel width as the 
characteristic length. The non-dimensional quantities come 
out to be, 


7 = Xg4, 

u(y)= 


(4.2.8) 
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ii(y) 


2y - y^- (l - cosh a y + tanh a sinh a y) 




(l “ l/cosh a) 


(4.2.9) 


The non-dimensional parameter K takes the form 
K = I * 2/^2 

= h2/£ 2 = ^2 (4.2.10) 

The governing equation for the stability of plane Poiseuille 
flow, sub;3ect to two dimensional disturbances, is obtained 
from equations (3.5.9) by substituting for values of K and 
u as, 


(d2-c>c2)2 


(d2- 
“ ^2 






0 = iC(R .S(u-c) (D^- £K^) 0 


(D^ u) 0‘ 


The boundary conditions are, 

0 = 0 

D0 = 0 at y = 0 and y = 2 

DS = 0 

where S = (D^- r/.^) 0 


(4.2.11) 


(4.2.12) 


Our task is to obtain a numerical solution of 
equation (4.2.11) with the boundary conditions (4.2.12), Vie 
are mainly interested in finding out the eigenvalues c and 
eigenfunctions 0 for a particular set of parameters R, (X.and a- 
A discussion of the general nature of the numerical solution of 
equation (4.2.11) is given in the next section. 
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4,3 NUI'^ERICAL solution OF THE DIFFERENTIAL EQUATION; 


The disturbance equation with couple stresses is 
given by equation (4.2.11). The parameter a comes in due to 
the effect of couple stresses. The case of a tending to 
infinity corresponds to the classical problem, and the resul- 
ting equation is the well known 0 r r- Somme rf eld equation ; 


(D^ 0 = io<R (u - c) (D^- x^)0-(D^ u)0 


Subject to the boundary conditions, 
0 = 0 
D0 = 0 


at 7 = 0 and y = 2 . 


(4.3.1) 


(4.3.2) 


The numerical solution of the Orr-Sommerf eld equation has 
been given by many authors. A brief description of the 
numerical solution of equation (4, 3. l) will help in under- 
standing the nature of equation (4.2.11). Hence, we will 
start our discussion with a review of tbe numerical solution 
of the Orr-Soramerf eld equation. ■ 

In solving the Orr-Soramerf eld equation, one is 
faced with the task of solving a linear ordinary differen- 
tial equation of the form 

Lg 0 + R"^ L^ 0 = 0 (4.3.3) 

where L^ is a fourth order operator and Lg is a second order 
operator. If R were of order unity there would have been no 
real problem. However, in the present problem R is of the 
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order of 10^, and hence the equation is singular in nature. 
Sinde u (y) is symmetrical about the centre line^ it is 
advigahle to consider the solution irdm the centre to one 
of the Walls* The boundary Conditions can then be specified 
at the centre and at the Wall. There will be two linearly 
independent solutions of the equation satisfying the two 
central boundary conditions. An appropriate linear combina- 
tion of these solutions must be obtained to satisfy the two 
wall conditions. It is known from the asymptotic theory 
(tinl) that one of the solutions grows rapidly away from the 
centre, a numerical solution of this function (Reynolds ) 
for R - lOOOO (corresponding to Thomas’s tabulated results) 
showed a growth of the order of lO^from the centre to the 
wall. It is clear that a numerical generation of the second 
solution, which is linearly independent, will be difficult . 
Any slight round-off error will in effect throw in some 
small multiple of the growing solution, which is likely to 
dominate the second solution by the time the wall is reached. 
This Was indeed observed in an experiment (Reynold^®) in 
which an eigenfunction calculation was started at the centre 
with the inviscid starting conditions. By the time the wall 
was reached the solution was, to eight digits, a multiple of 
the growing solution. In fact, it was of the order of 10^®, 
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suggesting that it was present initially to the order of one 

O 

part in 10 , the maximum accuracy of single precision on the 
machine. Of course, the final solution does not exhibit this 
rapid growth, which means that only a very small amount of the 
growing solution is required. The problem then is to generate 
numerically a second solution which is not merely a multiple 
of the growing solution. There are essentially two approaches. 
First, one might use multiple precision, extending the 
accuracy of the digital computations to, say, 16 digits. The 
solution can then be carried out in two parts, inwards from 
the centre and outwards from the wall and matched somewhere 
in between. This may limit the growth of the error to one 

O 

part in 10 on either side, for which 16 digit computations ^ 

might suffice. A scheme of this type was used by Nachtsheim. 

A second scheme used by Kaplan (1964) involves a suppression 

of the growing solution during the calculation of the well 

behaved solution. A general description of this technique can 

3 8 

be found in Reynolds paper. Kaplan' s method has been 
used in a variety of recent calculations with remarkable 
success. Experience has shown that the initial guess must 
be relatively good for fast convergence. 
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The finite difference technique suggested by Thomas 
and used by Kaplan is very time-consuming and complicated. 
Moreover, to start the numerical integration, the initial 
values. at considerably more points should bo known. Keeping 
these points in view, the method suggested by Kachtsheim has 
been used in solving eauatj.on (4.2.11). The technique will 
be discussed in detail in the following sectj.ons. 


4.40 METHOD OF SOLUTION ; INITIAL VALUE TEGflNlQUS 
The disturbance equation (4.2.11) is 


(D^ 


sf, 1 


1 - 




(D^ -0^2) 0 = I’ (u-c) 


(D^- (J?) 0 - (D^ u) 0' 


with the boundary conditions, 

0 = 0 

D0 = 0 at y = 0 and y = 2. 

DS = 0 


J 


The amplitude of the disturbance vorticity function is 
defined as . ‘ 

S = (D^ -o<^) 0 

If we introduce 

A = S 

then, equation (4.2.11) takes the form, 


A’ 


^ (d2-«2) ^ 

ar 


= i C>< R 



c) 


S-D^ 
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Thus, the sixth-order differential equation can be reduced 
to a system of- three second order differential equations, 

0 " = 0 + S ' 

S” = c<^ S + A '(4.4.1) 

A*’ = (a^ f’ - itX R (u - c) S - u "<56 

The corresponding boundary conditions are, 

0 = 0 

0' = 0 at 7=0 and y = 2 (4.4.2) 

S’ = 0 

where ' denotes dif f erent j.ation with respect to y. 

The solution of the disturbance equation (4.4.1) for 
0, for given values of ■>< , R and a, can be made to satisfy 
the boundary conditions (4.4.2) only for the ej.genvalues of 
c. We are interested in finding the eigenvalues and eigen- 
solutions for a given set of flow parameters. It is also 
of interest to determine the minimum critical Reynolds 
number. 

With regard to equation (4.4.1), of primary 
Interest are the solutions that are even functions of y about 
the line y = 1. Since the velocity profile is an even 
function of y about y = 1 , the disturbances can be separated 
into even and odd parts. The former, which has a simple 
flow pattern, usually give a lower critical Re^molds number, 
hence the second set of boundary conditions (at y = 2 ) are 
replaced by the conditions at centre, as 
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0 ' = 0 

0 at y = 1 (4.4.3) 

S* = 0 

Thus our characteristic value problem may be treated in the 
interval (0,1), in which u (y) is increasing monot eric ally. 

The basic idea is to solve the differential equation 
for assumed values of the eigenvalue and other initial values 
and to determine the appropriate changes in these values so 
that the boundary conditions are ultimately satisfied. 

Equation (4.4.1) basically has six linear independent 
solutions, some of which grow exponentially at a rapid rate. 
The difficulty is overcome by using Double-precision arith- 
metic (16 significant figures) and carrying out the integra- 
tion of equation (4.4.1) in two parts, inwards from the centre, 
y = 1, and outwards from the wall at y = 0 and then matching 
the two solutions in the middle at y^ = 5. 

In each case, the integration is performed in the 
direction in which the wanted function is increasing and the 
unwanted function is decreasing. It is important to include 
as much information as possible about the wanted solution in 
the statement of the problem. 
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In general the most important criteria are the proper 
choices of and Y (which in this case is Y = l). A poor 
choice of y^ could lead to convergence to a different solution. 
The best choice of y^ is that value of y for which, for the 
required eigenvalue, the coefficients of the dominating 
functions vanish in the equation which involves the derivatives 
of these dominating functions. If the chosen is somewhere 
in between the critical values corresponding to two eigenvalues, 
we may converge to either of these eigenvalues. Indications 
are that we would converge to the root closest to the initial 
guess but we may loose the desired speed of quadratic 
convergence. 

A poor choice of Y has two effects. The first is 
that if Y is too small to represent infinity for the eigen- 
solution sought, we may converge to a different result or 
may not converge at all. The second effect is that, if Y is 
much larger than necessary, the solution will converge but at 
a linear rather than at a quadratic rate. This is probably 
due to the fact that the best matching point is not the same 
for all functions and if we necessarily cover too large a 
range there is a significant build up of errors which 
inhibits the process. The shape of the eigenf •unction 0(y), 
but not the amplitude, is fixed by equation (4.2.11). There 
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are two approaches for fixing the amplitude of the eigen- 
function 0 (y) . 

Firstly, if we define the amplitude of the velocity 
distribution in some particular manner, this definition will 
fix the amplitude of 0, apart from a constant of modulus 
unity. An alternative approach is to normalise 0 in some 
manner. This will in turn define the amplitude of the velocity 
of the disturbances. In plane Poiseuille flow the latter 
choice is more convenient. Following the normalisation 
suggested by Nachtsheim we set, 

0(1) = 1 for even modes 

0‘(1)= 1 for odd modes 

As we are interested in even modes, we take the normalising 
condition 0(1) = 1. This fixes the size of the whole solution. 
Then our task is to solve equation (4.4.1) with the boundary 
conditions given by equation (4.4.2). As discussed earlier, 
the initial values at y = 1 and y = 0 are specified in the 
following manner. 

For the forward solution, starting from y = 1, the 


initial conditions are, 
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0f . — 1 

0 ^ = 0 

Sf ~ p (4.4,4) 

= 0 
Af = q 
A^ = 0 

The backward solution is started at y = 0 with the initial 
conditions, 

= 0 

0^ =0 (4.4.5) 

Sb = ^ 

s; = 0 

here the suffixes f and b represent the forward and the back- 
ward solution respectively. 

The variation parameters p and q in the forward 
solution and r, g and t in the backward solution cannot be 
fixed arbitrarily but must be determined in the iterative 
process whilst attempting to match the solutions at yc “ ‘S* 
The solutions must be continuous. Matching at y = y^ 
requires that 
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= 



= ■ ■ 


Sf 

= Sfc 


s' 

= sJ 

b 

( 4 ♦ 4 ♦ 6 ) 

Aj. 

^b 



= 



If these conditions are satisfied, all the higher derivatives 
agree and matching is accomplished. 

The quantities 0^ (y^) , 0 ^( 7 ^), SpCy^,), 8 ^( 7 ^), AfCy^.) 
and ApCyg) are functions of p, q and c and the quantities 
0 l 3 (yc)? -^b^yb^ functions 

of r, g, t and c. The quantities 0, S, A, p, q, r, g, t and 
c are complex. Successive changes are made in the first 
estimates of the variational parameters so that equations 
(4.4.6) are ultimately satisfied. 

The Newton-Haphson nethod is used to fulfil the 
conditions imposed by equations (4.4.6). If the chosen values 
of p, q, r, g, t and c produce a solution that approximately 
satisfies equation (4.4.6), a better approximation is obtained 
by starting with p+AP?l + 4iq, r+Ar, g + Ag, t + At 
and c + Ac. The increments Ap, Aq, Ar, Ag, At and Ac 
are solutions of the equations 



52 


0f - 0-b +Ap 

3 

3P 

+ Ar 

d 

3r 

0f - 0^ +Ap 


X U 

3P 

+ 

3 


3r 

Sf “ +AP 

3 

3p 

+ Ar 

3 


31* 

sl - Sv + Ap -A 

+ At 

a 


3r 

Af - A-b +Ap 

3 

SP 

+ Ar 

a 


ar 

^f ~ -^b 

3 

ap 

+ Ar 

a 


ar 


(0f - 01,) 

(0f - 0^) 

(0^ - 0^) 
(0f - 01,) 

tSf - Sb) 

(Sf - s^) 

- S'J,) 

(A^ — A-j^) 

(Af - Aj,) 

(Af - 4 ) 
(A- - aP 


+ A q C0f " 0b) + ^ (0f " 0b) 

+ A g (0f - 0|)) + At (0. - 0. ) = 0 

+ A q — ^ (0f - 0v) + AC (01 - 01 ) 

+ A g (0^ - 0^) +At ^ (0^ - 0^) = 0 

* ^ <j - Sb) +AC ^ (Sj. - Sb) 

+ A g ^ CSf - Sb) + At -A (Sf - Sb) = 0 

+ A<I^(S^ -Si)^Ae-A(s' -S;) 

" ■ ®b^ - S^,) = 0 

+ Aq ^ (Af - Ab) +AC ^ (Af - Ab) 

- Ag ^ (Af - Ab) +At (Af - Ab) = 0 

'4 ■ 4^ ^4 ■ 4^ 

+ Ag _A (a^ - Ap + At -A (Aj - aI) = 0 
"^S ot 

(4.4.7) 


in which the functions and the partial derivatives that constitute 
the coefficients are evaluated at y„. 

The partial derivatives are obtained by solving 
initial-value problems. These equations are obtained by 
partial differentiation of the terms in equations (4.4.1). 
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The coef f icj.ent s of equations (4.4.1) are analytic functions 
of y and the parameters , ’i, a and c. The solutions of 
equation (4.4.1) therefore have the same anal.vt.ic properties 


and posses the required partial derivati.ves . 

The quantities --5-™ = 0^ , = Sf- „ and — -~ 

f?P ^P 

= for the for.Arard solution satisfy the svstem of eouations 

1 5P 


0p 




5P 


+ s 


f sP 


S 


it 

f ,p 




(4.4.8) 


it 




+ a^) Af^p 


ic< Ra^ 



0 ) Rj.^p-u"0f^p 


with the initial conditions at y = 1 given by, 


f ,p 

^f,p 

,p 

^f,P 


^f,P 


A 


f ,p 


0 

0 

1 

0 

0 

0 


(4.4.S) 


Similar equati.ons can he written for the var.lational parameter 
q by replacing p by q in equations (4.4.8), and by satisfying . 
the following initial conditions at y = 1 " 
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%,q ~ ^ 




fvq 

t 

f ,q 
■f ,q 

t 


= 0 
= 0 
= 0 
= 1 
= 0 


(4.4.10) 


The variational equations for c are given by, 


%,c 

,0- ^ 

%,c 


q" 

^f,c 

y,c " 

■‘‘^f , c 


I 5 o 

= (a^ +o<.^) 

^ - io<R 

I ^ c 

(u 


(4.4.11) 


- c)5. 


u 0 


f , c 


With initial conditions at y = 1 given by, 


^f,c 


'f ,c 


= 0 
= 0 


%,c == 0 




= 0 
= 0 


f ,c 


( 4 , 4 . 12 ) 


For the baclo.'^ard solution, the variational equations for r, g 
and t are the same as those for p in the foward solution, 
except that the j.nitial conditions at y = 0 in this case are 
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(4.4.13) 


(4.4.14) 


for t and c the initial conditions at y = 0 are 


(4.4.15) 


1 
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For c, the initial conditions at y = O are 

c “ 0 

0 

Sh,c = 0 

Si’ = 0 

( 4 . 4 . 16 ) 


The quantities In the forward solution are functions of y, 

P, q and c only. Similarly In the backward solutions the’ 
quantities are dependent on y, r, g, t and o only. Hence 
the partial derlyatlves with respect to r, g and t in the 

forward solution and with respect to p and q In the baclrward 
solution are zero . 

That is, 




%,t = 

4,r 

= 0. 

D, q 

= 0 

,r 

= = 

^f,t = 



= 0 


= 4, a ■= 

•^f,t = 

'^b.p 

" '-1,0 

= 0 

4,r 

= o' 

f^g " 

q' 

^f,t ~ 


" ^b,q 

= 0 



h,t ^ 

^h>P 


= 0 

Ai 

‘ n , r 

= A i = 

f?g 

A ' 

,t 

4,p 

- 

= 0 
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Equa.tlons (4.4.7) then reduce to the following set of 
equations. 


0 f 

- S^b 

+ AP 

%,p 

q 

^f,q 

+ Ac 

«f,o 

- ^,c> - 

Ar f*b,r 







-AS 0 t,g 

“AP ^b,t 

= 0 

^f 

- 0’ 
^b 

+ Zip 

^f,P 

+Aq. 

0 l 

15 q 

+ Ac 


- 0b, c) - 

h',r 







- As 

^b,g 

0i,t 

= 0 

Sf 

- Sb 

+Ap 

Sf,P 

+A q 

%,q 

+ A c 

^Sf,c 

~ 5 b,C^ - 

A^ '^b,r 







-AS 

Sb,g 

"AP Sb^b 

(4.4.18) 

= 0 

S:^ 

- -^b 

+ Ap 

Sf,p 

+Aq 

4 ,q 

+ A c 

^®f,= 

- - 

Ar 3’^^ 







- Ag 

s,! „ - 
b,g 


= 0 

Af 

" ^b 

+ Ap 

f ,p 


%,q 

+ Ac 


■^b,c^ 

Ar "^b,r 







"AS 


-At Ab,t 

= 0 

4 

- ^b 

+ A p 

■^f ,p 

+ A q 


+A c 

‘Af,e 

- - 

Ar ' 







-Ag 

K,g 

-At 

= 0 


Hence, there are 6 complex equations to determine the 6 
complex quantities Ap, Aq, ac, Ar, Ag and At. 

Each step of the integration scheme is carried out 
hy starting with the estimates for p, q, c, r, g and t. The 
forward system of equations (4.4.1) is then integrated, step 
by step, sub.iect to the Initial conditions given by 
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equation (4.4.4). This will, necessitate the simultaneous 
solution of the variational equations (4.4.8) and simi.lar 
equations for q and c together with the in.ltial condj.tions 
given by equations (4.4.9), (4.4.10) and (4.4.12). The 
bacta/j'ard system is then integrated wj.th the initial condi- 
tions given by equations (4,4.5). The four variational 
system of equations similar to eauations (4.4.8) and (4.4.11) 
but with the initial conditions given bv equations (4.4.13), 
(4.4.14), (4.4.15) and (4,4.16) are integrated simultaneously 
with the backi^rard system of equations similar to equations 
(4.4.1). The forvrard and backward solutions are compared at 
the matching point, and the coefficients in equations(4.4 .18) 
are evaluated. Equations (4.4.18) are then solved forAp, 

A q, Ac, Ar, Ag and At, and the resulting solution gives an 
estimate for the Increments required for the next solution. 
This procedure is repeated till two consecutive solutions 
agree to a specified accuracy. 

Only variations with respect to the real parts of p, 
q, c, r, g and t.need to be obtained by this procedure. The 
solutions of equations (4.4.1) are analytic functions of 
p, q, c, r, g and t. Therefore, the real and imaginary 
parts of the complex derivatives, appearing in the coeffi- 
cients of equations (4.4.18), can be expressed in terms of 
the derivatives with respect to the real quantities only. 
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4.5 EQUATIONS IN REAL EOM : 


The differential equations, for the real quantities, 
are obtained by separating the original equations into real 
and imaginary parts. The equations (4.4.1), written in real 
form, come out to be, 


0^ = 0^ + 

0j = o(^ 0^ + 

S” = Sj. + Ay 

2 

S- = rV S. + A. 

1 ^■^11 

Tt / P P O 

A = (a + (X •) A + |5(R 

r r 


^i 


(U - c^) 


axid 


■>> _ 


(a2 + <j^2) g2 


- u 0. 


-S^ - Sy (U - Cy) 


if ^ 

+ u 0, 


(4.5,1) 


where the suffixes r and i represent the real and imaginary 
parts respectively. 

The real form of the variational equations (4.4.8) and 
and (4.4.11) comes out to be : 
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fr.Pr 




0f> — /V ^ 0J^ 

^^i,Pr ^^l,Pr ^^^ijPr 

®^r,Pr " '^^r,Pr %,Pr 


r,Pr 


s;i 

fi .P 


i,>^r 


^f.r. + o(^ S. 

^i,Pr ' %,Pr 


(4.5.2) 


''fr,Pr = ^ p^ ^ 


+ Sf. -n “ ■^'’ ^ 

J-i,Pr ^ 


^i,Pr 


and = (a2 + ot") c, 

- (u - =P u" 

The real form of the variational equations for q, r, g and t 
can he obtained from equations (4.5,2) by replacing p by q, 
r, g and t respectively. The variational equations for c, 
for both the forward and backward solutions, come out to be; 

= Sr,cr " *^r,C3, 


'"i c ~ ®i c o 

S" = A + j/ ^ p 

^ = A. + OC^ s. 

— Ti-vSi^/SN a jL_ 


P 9 Cp 


' r ' r P 

(a^ +of^) A^^q^ +p(R [- £ 

+ Sa p (u - Cy.) - S^• - u" 0. 


r,Cj, 1 


^i,Cr 




^i,cr " \,Oy, ^i,Cj. °i 

- Sr.Oj, - °r) + Sr u" »r.cl 
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The real part of the linear equations (4.4.18) comes out to be, 


( 0 |£‘) . 


(0^) + 

D ,r 

(0. ) 
f ,p r 

APr - 



-- 


^i A qp 

*i4,c 



- 


oh - 



+ 


i - 


P A gj, + ^^b,g^iASj_ 

- 



^ ^^b,t 

). Atj = 0 



r Aq, 


Similarly, the other equations can also be split into their 
real and imaginary parts. 

Since the derivatives, with respect to the real 
quantities, are the one that are actually calculated, 3.t 3. s 
necessary to express the real and Imaginary parts of the 
complex derivatives in terms of the real quantities only. 
For example, in the case of 0f pi (0p^p)j. = ^^f^r,p 


Thus we have to solve twelve equations for the 
corrections and thirty variational equations for obtaining 
the desired quantities. 


4.6 INTEGRATION FORMULAS ; 

The integration is performed by using the fifth- 
order Predictor Corrector method of Milne, which uses a 
fourth-order Runge-Kutta method to obtain starting values. 
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Let the system of n equations, to be solved, be given 
in the form 


y^ = (x, y^, Yg 


y^^) (i = 1, 2 


.n) 


(4.6.1) 

with the initial conditions 

Ji (Xq) = y^^ •, yi (Xq) = y^^ (i = 1 , 2 ....... n) 

Let and be the values of and yj^ at x = Xj^ and 

^i,k second derivative of yj^ at x = X|^. h denotes the 

step size. The special Runge-Kutta formulas used are as 
follows s 

%! = “i 

Rs “ 2> n,k i n,k s 

h 


K. 


13 = + J'’ n.k * K,k i 5^13) 


( 4 . 6 .S) 


yi,k+r yi,k: 


h,k i 

12 ^is' 


and 


4,k.i= n,k"i 'hi" 


''here fi stands for f^ (Xj^, 

The Milne Predictor-Corrector formulas for solving the 
system, given by equation (4.6.1), are 

Pl,k+1 = n,k " U,k-2 - yi,k-3 " hy4 (Sfi,i,+2fi'^jj.^-Sfi^^_3) 

n,k+l “lyijk ■ n,k-l " hlis Pi^jj-.^) + lOU^jj: 




(4.6.3) 
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The Corrector formula, given by equation (4.6.3), is applied 
only once, so that only two derivative evaluations are needed 
for each step of the Milne integration. The starting values 
needed in the Predictor formula (4.6.3) are obtained by 
using equations (4.6.2). 



5 . RESULTS AND CONCLUSIONS 


5.1 MERIC AL RESULTS FOR PLANE POISEUILLE FLOW? 

To study the stability of plane Poiseuille flow^ the 
eigenvalue c have been calculated for different flow parameters 
Relations between various flow parameters have been shown 
graphically. 

Couple stresses affect the- steady state velocity 
profile. It can be seen from Fig. 1 that at any point, the 
steady state value of the velocity is smaller for lower values 
of a. This essentially means that the effect of couple stress 
is equivalent to an apparent increase in viscosity. 

Neutral stability curves for plane Poiseuille flow for 
different a have been shown in Fig. 2. Part of the neutral 
stability curve (Lin^) for the nonpolar case has also been 
plotted for purposes of comparison. The neutral stability 
curve Cj_ (»<, R, a) - 0. separates the regions of stability 
and instability in the -R plane. The lower branch of the 
stability curve exhibits the stabilizing nature of viscosity. 

On the upper branch, the effect- of viscosity is to destabilize 
the flow. The effect of this destabilizing mechanism is to 
shift the phase difference between the u' and v’ components of 
the disturbances. This helps in the development of the 
Reynolds stress - ^u'vV, proportional to (010^ - 0r0i)5 which 
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is responsible for the enefgy transfer between the meah flow 

and the disturbances; If at any stage, the energy supplied by 

» 

the Reynolds stress becomes greater than the energy dissipated 
by the viscosity through the diffusion of vorticity^ the dis- 
turbances are amplified; For smaller Reynolds numbers the 
effect of viscosity is to stabilize the flow. The minimum 
critical Reynolds number indicates the value of the Reynolds 
number below which the disturbances are damped. It can be seen 
from Figs. 2 and 3 that the minimum critical Reynolds number 
first decreases and then increases with the increase in the 
couple stresses. Also from Figs. 2 and 4 it follows that the 
critical wave number increases with the increase of couple 
stresses (that is with decreasing a).' The limiting case of 
zero couple stresses (a— corresponds to the equation of 
motion (3.5.3) which governs the stability in the nonpolar case. 
It can therefore be assumed that the stability curve for a — >00 
will be the same as the stability curve for the nonpolar case. 
From Fig. 2 it can be seen that, in the presence of couple 
stresses' the domain of the unstable region increases appreciably 
in comparison to the nonpolar case. It should also be noted 
that the instability occurs at higher wave numbers in the 
presence of couple stresses. Table 1 gives the minimum criti- 
cal Reynolds numbers and the corresponding wave numbers for 
different values of a as well as for the nonpolar case. 
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The mode shapes of the eigensolutions 0 for the wave 
number = 1 and Reynolds number R=6000j have been shown for 
different values of a in Fig, 5. It can be seen that 0^ is on 
the whole much greater than 0 ±, It can also be seen that 0j_ 
is concentrated in the critical layer region (where the phase • 
velocity Cp equals the mean velocity u(y))near the wall. The 
Reynolds stress and the couple stresses are also found to be 
concentrated in the critical layer region. When the distur- 
bances are damped,, they will vary steeply in the critical 
layer. This can be seen from Fig. 5 in which the mode shapes 
for a = 5 and a = 10 correspond to damped oscillations whereas 

those for a = 15 correspond to amplified oscillations. The 
rapid oscillations in the critical layer regions are truly 
remarkable and give rise to the diffl culty in direct numerical 
integration. 

Figs. 6, 7 and 8 show the variations of the eigenvalue 
Ci with the wave number c< for different values of R and a. For 
sufficiently low values of R the flow remains stable. Corres- 
ponding to each Reynolds number, there are two critical wave 
numbers for which Cj[=0. These two values define the range of 
instability. 

The variation of Cj_ with R for different values of 
and a can be seen from Figs. 9, 10, 11, 12 and 13, For any 
fixed , the critical Reynolds number corresponds to zero 
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value of c^. Figs. 3 and 4 show the variation of the -minimuin 
critical Reynolds number and the minimum critical wavd number 
with a. 

5.2 KE2^ ARKS ON THE NUMERICAL SOLUTION; 

The instability is predicted for large values of the 
parameter cARa^. The per step truncation error involved in the 
numerical integration increases with o^Ra^. To have a better 
control on the per step truncation error an appropriate step 
size must he chosen. This was achieved by trial. With a 
step size of 0,0025 two iterations can be performed in one 
minute. It is essential to have good approximations for the 
initial values for fast convergence. Tables 2, 3 and 4 give 
the initial values at representative points for different 
values of a. Lagrangian interpolation was used for determining 
the critical Reynolds numbers and the critical wave numbers by 
using the values obtained for different values of the flow 
parameters. 

Appendix B contains a list of symbols. A brief dis- 
cussion of the program is given in appendix C. A listing of 
the program has been given in appendix D. 

5.3 CONCLUDING REMARKS: 

The couple stress theory introduces a material 
constant ^which has the dimension of length. A size dependent 
effect is predicted by this theory which is not present in the 
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nonpolar case. The effects of couple stresses are quite large 
for small values of the parameter a = h/f • The parameter a 
enters the resulting equations of motion as the multiplier of 
the highest order space derivative. The limit a , 
corresponds to the classical nonpolar case. The order of these 
equations and the number of boundary conditions is diminished 
in the transition to the classical theory. 

From the discussions in the previous sections it follows 
that the effect of couple stresses is equivalent to an apparent 
increase in the viscosity of the fluid. 














FIG 6 - EIGENVALUES FOR PLANE POISEUILLE FLOW FOR 
VARIOUS R. 
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FIG. 7 -EIGENVALUES FOR PLANE POISEUILLE FLOW 
FOR VARIOUS R. 




FIG. 8- EIGENVALUES FOR PLANE POISEUILLE FLOW 
FOR VARIOUS R. 
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FIG.!3- EIGENVALUES FOR PLANE POISEUILLE FLOW 
FOR VARIOUS . 
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APPENDIX B 


A 

a 
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h 
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t 

u 

v’ 
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c< 
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Symbols 

S"- c<2s 

non-dimensional couple stress parameter 

phase velocity 

A(0) 

half -width of the channel 
curvature-twist rate tensor 
material constant 
S(l) 

A(l) 

S(0) 

Reynolds number 

disturbance vorticity amplitude 
A«(0) 

velocity of basic flow 

disturbance velocity parallel to plates 

disturbance velocity normal to plates 

distance parallel to plates 

normal distance from lower plate 

matching point 

wave number 

vorticity vector 

stream function amplitude 

stream function 



Superscripts ; 


ix 


* dimensional quantities 

’ denotes differentiation with respect to y 

(--) basic flow quantities 

(.) time derivative following a particle 



APPENDIX G 


DESCRIPTION OF THE FORTRAI-I PROGRAM FOR THE SOLUTION OF THE 
EIGENVALUE PROBLEM FOR PLANE POISEUILLE FLOW 

The ntmierical procedure outlined previously for 
solving the eigenvalue prohlem was programmed for solution on 
the IBM 7044 in fortran IV. 

The correspondance between the FORTRAN symbols used 
in the program and the mathematical notation employed previou- 
sly is shown in the following list ; 


FORTRAN^ 

symbol 

Mathematical 

symbol 

fortran 

symbol 

Mathematical ' 
symbol 

Y1 

^2*r 

YIC 

^r,Cp 

Y2 

0i 

y2C 

'*'1 c 

SI 


SIC 

^r c 

S2 

Si 

S2C 

®i c 

A1 

Ar 

AlC 

A 

^r,Cr 

A2 

Ai 

A2C 


YIP 


yiq 


Y2P 


Y2q 

^ 1 , 1 , 

SIP 

^?Pr 

1 SlO 


S2P 

%,p^ 1 

S2Q 

S j 

^jlr 

AlP 


AIQ 


A2P 

^ijPr 

A2Q 



xi 


FOEERM 

symbol 


Mathematical 
gymb ol 


®i,t 



FORTEM 

symbol 

DAIG 

DA2C 

DYIQ 

DY2q 

DSIQ 

DS2Q 

DAlQ 

DA2Q 

DYlT 

DY2T 


DSlT 


Mathematical 

svmbol 


Tj Cj. 


^i<lr 

0 ' 

3-jqj. 

S’ 

r,qr 
Si ^ 

^ ? q^ 
K ^ 

^jq^ 

Al 

ijqp 

t 

0 ’ 

i.t 

7 “•»» 


DYIP 

DY2P 

DSIP 

DS2P 

DAIP 

DA2P 

DYIC 

DY2C 

DSIG 

DS2G 


0' 

r,Pr 

01 

S’ 

TjP j. 

^i?Pr 

^^jPr 

^i?Pr 




DS2T 

BAIT 

DA2T 


DELPl 

DELP2 


DEL ni 


DELQ2 

DELRl 

DELR2 


SI 

ijtj. 

^r,ty 

H. t 


z^qr 




xii 


FORTRAN 

SYmlDOl 

Mathematical 

symbol 

fo^ran 

symbol 

Mathematical 

symbol 

DELGl 


aiback 

ApCO) 

DELG2 


A2BACK 

A^CO) 

DSLTl 

^tj. 

daibak 

A'(0) 

CELT 2 

4ti 

DA2BAK 

Ai(0) 

DELCl 


AL 

c< 

DBLC2 

A 

R 

R 

SIFTO 

Syd) 

AHL 

a 

S2FWD 

3^(1) 

W 

u 

AlFWD 

Ay(l) 

DDW ' 

u” 

A2F¥D 

A^(l) 

X 

7 

sibagk 

1 SyCo) 

XSND 

7=1 

S2BACK 

Si(0) 

XMATCH 

7c 


xili 

The following remarks are intended to aid in a study 
the program t 

(1) Subroutine JEJ is used to evaluate the second 
derivatives. The variables 2 and DDZ -that appear 
in JEJ are dummy variables. 

(2) Subroutine JAlh is used to store the matrix of 
coefficients that are formed from functions and 
partial derivatives evaluated at the matching 
point. The solution of t-he simultaneous linear 
equations is accomplished by calling subroutine 
Gauss. Before the subroutine GAUSS is called, EE 
contains the coefficient matrix and W contains 
the ^’right-hand side”. After GAUSS is called, W 
contains the answers. 

(3) Subroutine lUTEGR carries out the step-by-step 
integration with either the Runge-kutta method 
(INDEX=0) or the Milne n^hod, which uses the Bunge 
kutta method to obtain starting values (IUDlIX=l). 
The program listing is given in appendix D* 



o o n o o 


APPENDIX D 


£JOB ME6092,TIME008,PAGES030,NAME JKJAIN 
ilBJOB 

£IBFTC JKDP NOPRNT 

STABILITY OF PLANE POISEUILLE FLOW 


EXTERNAL JKJ 

DIMENSION T(30),DTC30J,DDT{30)tSS(200)tVK{l5) »DELC16I ,VVKC16> 

DIMENSION 0eLT(L5,15O) 

DOUBLE PRECISION 

1T,DT, SlFW0,S2FWD#AlFWD,A2FW0,SiBACK,S2BACK,AlBACK,A28ACK,DAlBAK» 402 
1DA2BAK,DELPI,0ELP2,DELQ1,DELQ2,DELR1,DELR2,DELG1,0ELG2,DELT1,DELT2 403 
2,DELC1,DELC2*C1,C2 

DOUBLE PRECISION Yl, Y2, SI, S2, Al, A2, YIP, Y2P, SIP,S2P» AlP, A2P, 
lYlC,Y2C,SlC,S2C,AlC,A2C,YlQ,Y2Q,SiQ,S2Q,AlQ,A2Q,YlT,Y2T,SlT,S2T, 
2A1T,A2T,DY1,0Y2,DS1,DS2, DAI, DA2,0Y1P,DY2P, OSIP, DS2P,0A1P,DA2P, 
3DY1C,DY2C,DS1C,DS2C,DA1C,DA2C,0Y1Q,0Y2Q,0SIQ,DS2Q,DA1Q,DA2Q, 

4DY IT , DY2T , DS IT , DS2T , D A IT , DA2T 
DOUBLE PRECISION DDT, ACY1,ACY2,DEL,VK,VVK,SS,SFRESH,DELT 
COMMON T ,DT 

COMMON S1FWD,S2FWD,A1FWD,A2FWD,SIBACK,S28ACK,A1BACK,A2BACK,0AIBAK, 402 
1DA28AK,DELP1,DELP2,DELQ1,0ELQ2,DELR1,DELR2,0ELG1,DELG2,DELT1,DELT2 403 
2,DELC1,DELC2 

COMMON C1,C2,AL,R,W,D0W,AA,AR,AHL 401 

EQUIVALENCE ( Y1,T 1 1 T) , ( Y2,T{ 2 J I , C Sl,T (3) ) , I S2,TI4 ) I , ( A1 ,T ( 5) I , ( A2 , 406 
IT (6) ) ,(YlP,T(7)),{Y2P,Tl8n, (S1P,TI9) ),(S2P,TC 10) ) , { A1P,T(11) ) , { A2 407 
2P,T{12) ) ,(Y1C,T( 13) ),{Y2C,T{14)),(S1C,T{ 15) ),(S2C,T(i6) ),(A1C,T(17 408 
3) ),{ A2C,T{18) ),(Y1Q,T( 19 ) ) , t Y2Q,T( 20) ) , ( SlQ,Tt 21) ) , I S2Q,T ( 22) ) , ( A1 
4Q,T(23)) ,( A2Q,Tl24)),(YiT,T{25)),tY2T,T(26J),ISlT,T(27) ),<S2T,T128)), 

5) ),( A1T,T(29) ) ,( A2T,T(30) ) 

EQUIVALENCE { DYl , DT ( 1 ) ) , { DY2,DT{ 2 ) ) , C DSl , DTI 3 ) ) , (DS2, 0T{4) ) , 

1(DA1,DT(5) ),(0A2,0T(6)),IDY1P,0TI7) ) , IDY2P,DT(8 ) ) , (DSIP,DT(9) ) , 
2(0S2P,DT(10)), (DA1P,DTI II) ) , { DA2P,DT( 12) ) , IDY1C,DT( 13 ) ) , IDY2C,DTI 14) ), 

34) ), IDSlCtDTI 15) ), (DS2C,DT(16)), {DA1C,DTI 17) ),tDA2C,DTfl8) ), 

4(0Y1Q,DTI 19) ), (DY2Q,0TC20) ),(DS1Q,DTI21)) ,(DS2Q,DT(22)I, 

5(DA1Q,DT(23)), (0A2Q,,DT{24)),(DY1T,DT(25)) ,(DY2T,DTI26)),{DSIT,DTI27) ),I OS 
67) ), I DS2T,DT{28) ) ,(DA1T,0T(29) ), {DA2T,DT(30) ) 

C *******»«.*****«.****.*****»*♦*»********»*»»♦***■»*****»*******»*♦#***»»•**««■* 

READ 70, INOEX,NFWD,NBACK,NR,NAL,NAHL,ITERAT 
70 FORMAT! 715) 

READ 71,H,DELX,XEND,XMATCH 420 

71 F0RMATI4F6.4) 

READ 72,ACYl,FR,FTN,IPRINT 

72 F0RMAT(3F8.4, 15) 

READ 73,HR,HAHL,HAL 

73 F0RMATt3F14.6) 

READ 74,AL,AHL,R 
FQRMAT(F12.8, F12.8,F12.4) 


74 


425 


READ75,SlFWD,S2F»iD,AlF»iD,A2FK0,Sl8ACK,S2BACK,AlBACK,A2BACKfDAl8AK, 428 
1DA2BAK,C1,C2 429 

75 F0RMAT(3D24.16I 
PRINT 130 

PRINT 76» INDEX, NFWD,N8ACK, NR, NAL,NAHLtITERAT 

76 F0RMAT{//2X,*IN0EX=*,I5,2X,*NFWD=*, I5,2X,*NBACK=*,I5,2X,*NR=*,I5, 

1 2X,*NAL=*, I5,2X,*NAHL=*, I5,2X,*rTERAT=», 15) 

PRINT 77,H,DELX,xeN0,XMATCH 434 

77 F0RMAT(//2X,3H H=Fl2.8,2X,5H0eLX=F12.8,2X,5HXEND=F12. 8,2X,7HXHATCH 435 

1=F12.8) 435B 

PRINT 78,ACY1,FR,FTN,IPRINT 

78 F0RMAT(//2X,*ACYI=*,F12.8,2X,*FR =»,F12. 8,2X,*FTN=», F12. 8, 

1 2X,*IPRINT=*, 15) 

PRINT 79,HR,HAHL,HAL 

79 F0RMAT(//2X,3HHR=F8.2,5X,5HHAHL=F12.8,5X,4HHAL=F12.4) 

PRINT 130 

C ************************************************************** ********•»**»* 


80 


81 

C 


c 

110 


49 


FTC=FTN/5.0 

FRC=FR 

FTNS=FTN 

IPRNT=IPRINT 

IRT=ITERAT 

CALL FLUN(5000) 

CALL FLOV(5000) 

AST=0, 440 

BST=0. 441 

CST=0. 


00 lOO KAHL=1,NAHL 

AHL=AHL+BST 

DO 100 KAL-1,NAL 

AL=AL+CST 

DO 100 KR=l,NR 

R=R+AST 

ITERAT=IRT 445 

PRINT 80,AL,AHL,R 446 

F0RMAT(/2X,3HAL=FI2-8,2X,4HAHL=F12.8,2X,2HR=F12.4/) 447 

PRINT 8l,SlFWD,S2FW0,AlFWD,A2FWDtSl8ACK,S2BACK,A18ACK,A2BACK,DAlBAK448 
1K,0A2BAK,C1,C2 449 

FORMAT! 2X,3D24. 16) 


AA=AL-»*2 

AR=AL*R 

JK=l 

IFINAL=2 


454 

455 


1 = 1 
J=1 

JK=JK+1 

CONTINUE 

Y1P=0. 

Y2P=0. 

DY1P=0. 

DY2P=0. 

S1P=1. 

S2P=0. 

DS1P=0. 


456 

457 

458 

459 

460 

461 

462 

463 

464 

465 


DS2P=0. 

A1P=0. 

A2P=0. 

DA1P=0. 

DA2P=0. 

Y1C=0, 

Y2C=0. 

0YIC=0. 

DY2C=0. 

S1C=0. 

S2C=0. 

DS1C=0. 

DS2C=0. 

A1C=0. 

A2C=0. 

DAIC=0. 

DA2C=0. 

Y1Q=0. 

Y2Q=0- 

DY1Q=0. 

DY2Q=0. 

SiQ=0. 

S2Q=0. 

DS1Q=0. 

DS2Q=0. 

A1Q=1. 

A2Q=0. 

0A1Q=0. 

DA2Q=0. 

GO TO (50,51) tJ 
50 J=2 

C FORWARD SOLUTION 

N=NFWD 
M=1 

X=XEND 

H=-ABS(H) 

XPRINT=XENO 

DELX=-ABS(DELX) 

Yl=l.O 

Y2=0. 

DY1=0. 

DY2=0. 

S1=S1FWD 

S2=S2FWD 

DS1=0. 

DS2=0. 

A1=AIFWD 

A2=A2FWD 

DA1=0. 

DA2~0» 

GO TO 60 
51 J = 1 

C BACKWARD SOLUTION 

N=NBACK 
M=2 

Y1T=0- 


466 

467 

468 

469 

470 

471 

472 

477 

478 

479 

480 

481 

482 

483 

484 

485 

486 


t 

487 

488 

489 

490 

491 

525 

526 


528 

497 

498 

499 

500 

501 

502 

503 

504 

505 

506 

507 

508 

509 

510 

511 

512 

513 


PRINT 130 

130 F0RMAT(/2X, 120(lH-n 

C »***♦■»**»**«***»**»**»»»*»**«■*****»*»*♦«*■^^********»■***•»*»***•*»**•****^ 

DEL(1)=DELP1 
DEL(2)=DELP2 
DEL{3)=DELai 
DEL(4)=0ELQ2 
DEL(5)=DELR1 
DEL(6)=DELR2 
DEL(7)=0ELG1 
DEL(8)=DELG2 
DEL(9)=DeLTl 
0EL{ 10)=DELT2 
DEL111)=DELC1 
DELC 12)=0ELC2 
DO 1005 11=1,12 
1005 DELT( II, JK)=DEL(IU 

IFdFINAL.EQ.DGO TO 105 
SS(JK)=Y1 

IF( JK.EQ.2)G0 TO 475 
IF (SS{JK)-SS{JK-1 11475,375, 375 
475 VK(1)=S1FWD 
VK(2)=S2FWD 
VK(3)=A1FWD 
VK(4)=A2FWD 
VK(5)=S1BACK 
VK(6I=S2BACK 
VK(7)=A1BACK 
VK(8)=A2BACK 
VK(9)=DA1BAK 
VKl 10)=DA28AK 
VK(11)=C1 
VK(12)=C2 
DO 477 1=1,12 
477 VVK(n=VK(I) 

IPRINT=IPRNT 

JN=1 

ICOUNT=0 

FR=FRC 

FTN=FTNS 

IF(Y1-ACY1)105,105,275 

375 SS{ JK)=SSt JK-1) 

IPRINT=IPRINT+1 

JN=2 

IC0UNT=IC0UNT+1 

FTN=FTC 

IF{ IC0UNT.6e.2}FTN=FR 
IFl IC0UNT.GT.2)FR=FR/2. 

00 376 1=1,12 

DELTI I,JK)=DeLT( I, JK-1) 

OEL(I)=DELT(I,JK) 

376 VKin=VYK(I)+FTN»DEL{ 1) 

GO TO 405 

275 CONTINUE 

DO 276 1=1,12 

276 VK(I)=VVK{I) + FTN*OELtn 


405 CONTINUE 

S1FWD=VK{1) 

S2FWD=VK(2) 

A1FW0=VK{3) 

A2FW0=VK(4) 

SIBACK=VK(5) 

S2BACK=VK(6) 

A1BACK^VK(7) 

A2BACK=VKI8J 

0A1BAK=VK(9) 

DA2BAK=VK{ lOJ 
Cl=VK(ll) 

C2=VK(12) 

PRINT i25tITERAT,FTN, JN,FR 

125 F0RMAT(10X,*ITERAT =*,I5,20X,*FTN =»,F8.4,20X,*JN =*, I5,5X,*FR=*tF8.4 
I 8.4) 

PRINT 81,S1FWD,S2FWD, AlFWD,A2FWD,SlBACK,S28ACK,Al8ACKtA2BACK,DAlBAK448 
1K,DA28AK,C1,C2 449 

ITERAT=ITERAT+l 
IF(ITERAT-150)55,55,105 
55 IFtYl-ACYl)105,105,ll0 

105 AST=HR 581 

8ST=HAHL 
CST=HAL 
PRINT 140 

X40 F0RMAT{/2X,120{1H*)//20X,15HFINAL RESULTS//2X, 1201 IH*) ) 

PRINT 80, AL,AHL,R 

PRINT 8i,SlFWD,S2F«0,AlFWD,A2FKD,SlBACK,S28ACK,A18ACK,A2BACK,0Al8AK448 
1K,DA2BAK,CI,C2 
PRINT 54,Yl,Y2 
PRINT 130 
IFINAL=IFINAL-1 
IFdFINAL.EQ.DGO TO 110 

100 CONTINUE 583 

STOP 

END 585 

£IBFTC INTDP NOPRNT 

SUBROUTINE INTEGR(M,N,H,X, ISET,Y,0Y,D0Y, INDEX, JAI ) 2 

DIMENSION Y{30),0Y( 30), DDY( 30), YLLL 130), YLLC30) ,YL (30), YR(30) , 

10YR(30) , DYL(30),DDYLL{30) ,0DYL{30),0DYR{30) ,RK2(30),RK3{30) , 

2P(30) ,T{30),DT{30) 

DOUBLE PRECISION 

IT,0T, S1FWD,S2FWD,A1FWD,A2FWD,S1BACK,S28ACK,A1BACK,A2BACK,0A1BAK, 402 
lDA2BAK,DELPl,DELP2,DELQl,DELQ2,DeLRl,DELR2,0ELGl,DELG2,DELTl,DELT2 403 
2,DeLCl,DELC2,Cl,C2 

DOUBLE PRECISION E,YLLL,YLL,YL,Y,DYL,DY,ODYLL,ODYL,DDY,fR,DYR, 

1 D0YR,RK2,RK3,P 
COMMON T ,DT 

COMMON S1FWD,S2FWD,A1FWD,A2FWD,SIBACK,S2BACK,A18ACK,A2BACK,DA18AK, 402 
10A26AK,DELPl,DELP2,DELQl,DELQ2,OELRi,DELR2,DELGl,DELG2,DELTl,DELT2 403 
2,0ELC1,DELC2 

COMMON Cl, C2,AL,R,W,D0Sfli,AA,AR,AHL 2A 

E=H 5 

IF(ISET)30,30,31 
30 IF(INDEX)32,32,34 

32 K=1 


8 


33 

34 

31 

35 

38 

39 

41 

48 

60 

61 

62 

63 


42 

145 

36 

37 

43 

4 1 

45 


CALL JAI (M,X,Y,DDYJ 

GO TO 145 

K=2 

GO TO 33 

GO TO (35,36,36,36,37J,K 
DO 38 1 = 1, N 

P{ n = Y< n + ( H/2, }»DY{ I l+(H*H/8. )*DOY( I ) 

CALL JAI(M,X+H/2. tP,RK2> 

DO 39 1=1, N 

P{ n = Y(I )+H*DY{ I) + CH*H/2.)*RK2(n 
CALL JAI (M,X+H,P,RK3) 

DO 41 1=1, N 

YRtlJ =Y( n+H*COY( n + (E/6, J*{DOYC n+2.*RK2{ I J n 
DYR{ I )=DYC I ) + { E/6. )*{ ODYC I } +4.»RK2( I ) +RK3n ) ) 
CALL JAI{M,X+H,YR,D0YRI 
CONTINUE 
DO 42 1=1, N 

GO TO (60,60, 61, 62, 63), K 
Y{I) = YR(n 
DYII)=0YR{1) 

DDY(I)=DDYR( U 
GO TO 42 
YLL(U=0. 

YL(I)=0. 

ODYLC I)=0. 

GO TO 63 
YLL{n=0. 

GO TO 63 
YLLL(I)=YLL( n 
YLL(n=YL{I) 

YL( n=Y( I) 

Y(I)=YR(I) 

DYL{ n=DY{ n 

DY(I)=OYR(n 
DOYLL{n=ODYL( I) 

DDYL{I)=DDY(I) 


10 

11 

12 

13 

14 

15 

16 
17 

19 

20 
21 

23 

24 

26 


31 

34 


28 

29 


31 

32 

33 


DOYd )=DDYR(I) 

CONTINUE 
RETURN 
K=K+1 
GO TO 35 
DO 43 1=1, N 

p{ I ) = Y( I ) +YLL{ I)-YLLL ( I ) + ( H»H/4. )*( 5-*DDY ( I)+2.*DDYLt II-S.-^DDYLLC I M 

1) y 

CALL JAItM,X+H,P,DDYR) 

00 44 1=1, N 

YR{ I)=2.»Y( n-YLt n + (e*E/ 12 .)*(DDYR(I) + 10.*D0Y{ n+DDYLlI) ) 

CALL JAI{H,X+H,YR,ODYR) 

DO 45 1=1, N 

OYR( n=DYL(n+(E/3-)*(D0YR(I)+4-*0DY( I)+DDYL( n) 

GO TO 48 


35 

36 

37 

38 

39 

40 

41 

42 

43 

44 
45 

46 

47 

48 


END 

£IBFTC jkjdp noprnt 

SUBROUTINE JKJ ( M,X,Z, ODZ ) 

DIMENSION T{30>,DT{30),Z(30),DDZ(30) 
DOUBLE PRECISION 


IT ,DT, SiFWD,S2FWD,AlFWD,A2FWD,SlBACK,S2BACK,AlBACKf&2BACK,DAlBAK, 402 
10A2BAK,OeLPl,DELP2,OELQl,DELQ2,OELRl»DELR2,DELGl,OELG2,OELTl,DELT2 403 
2,DELClfDELC2,Cl,C2 
DOUBLE PRECISION Z,DOZ 
COMMON T ,-OT 

COMMON SlFWD,S2F«Dt AlFW0,A2FviD,Sl8ACK,S2BACK, AlBACK,A28ACK,DAlBAK, 402 
IDA2BAK,DELP1,0ELP2,DELQ1,DELQ2,D6LR1,DELR2,0ELGI,DELG2,DELTI,DELT2 403 


2,0ELC1,DELC2 

COMMON C1»C2,AL,R,W, DOW, AA»AR,AHL 102 

AHLX=AHL*X 104 

AAHL=AHL*AHL 105 

BP=AL*R*AAHL 106 

BNUM= 1--C0SH< AHLX )+TANH{ AHL )*SINH( AHlX) 

BDEN=l--2-/AAHL*C l.-l./COSH(AHL) ) 

BRATI0=BNUM/8DeN 113 

W=(2.*X-X**2-2./AAHL*BNUM)/BDEN 114 

DDW=-2.*BRATI0 115 

SUM=AAHL+AA 116 

0DZ{1I=AA*Z{ 1)+Z{3} IIT 

D0Z{2)=AA*Z(2>-»-Z(4} 118 

DDZ(3)=AA*Z{3J+ZC5) 119 

0DZ(4)=AA*Z{4)+Z{6} 120 

DDZ(5)=SUM*Z(51+BP*i-Z(3)»C2+Z(4l*(W-Cl)-0DW*Z{2y ) 121 

D0Z{6)=SUM*Z(6y+8P*{-Z{4)*C2-Z{3l*(«-Cil+DDW*Z(iy ) 122 

0DZ(7)=AA*Z(7)+Z(9) 123 

0DZ(8)=AA*Z(8}+Z( 101 124 

DDZ{91=AA*Z(91+Z( 11) 125 

ODZ(IO)=AA»Z(101+Z(12) 126 

DDZ{ll)=SUM»Z(ll)+BP*(-Z{9)*C2+ZC10l*(W-Cl)-DDW*Z(8)) 127 

ODZ(12)=SUM*Z(12)+BP*{-Z(10)*C2-Z(9)*(W-C1)+DOW*Z(7)) 128 

DOZl 13)=AA*Z( 13)+Z{ 15) 129 

0DZ(i4)=AA*Z{ 14)+Z(16) 130 

DOZ{ 15)=AA*Z( 15)+Z{17) 131 

D0Z{16)=AA*Z{16)+Z118) 132 

DDZ( 17)=SUM*Z( 17)+BP*(-Z( 15)»C2+Z{ 16)*{W-Cl)-DDW*Z(14)-Zt4) ) 133 

ODZ( ia)=SUM»Z ( 18)+BP*{-Z( 16)*C2-ZC15) *< W-C1)+0DW*Z{ 13 )4-Z(3) ) 134 

0DZ(19)=AA»Z( 19)+Z{21) 137 

D0Z120)=AA*Z(20)+Z(22) 138 

DDZ(21)=AA*Z(21)+Z(23) 

DDZ(22)=AA*Z(22)+Z{24) 1 

DDZt23) = SUM*Z{ 23)+BP«(-Z( 2l)*C2+Z(22)*(W-Cl)-0Dl!i*Zl2O)) 141 

DDZ(24)=SUM*Zl24)+8P*(-Z(22)*C2-Zt21)*(W-Cl)+D0W*Z(19)) 142 


GO TO t63,64),M 
6+ CONTINUE 

DDZ125)=AA*ZC25)+Z(27) 

DDZ(26)=AA*Z(26)+Z(28) 

DDZ(27)=AA*Z(27)+Z(29) 

D0Z(28)=AA*Z(28)+Z(30) 

DDZ (29) = SUM*Z( 29) +BP*( -Zt 27y*C2+Z(28)*( W-C1)-DDM*Z< 26) ) 
D0Z(30y = SUM»Z(30) +BP*l-Z{28)*C2-Z{27)*{ W-C1)+OOW*Z(25) ) 

63 RETURN 

END r 

£IBFTC JAINDP NOPRNT 

SUBROUTINE JAINt It IPRINT) 

DIMENSION T(30)»DT130),EE{ 12,12)fVV(12f 1) 

DOUBLE PRECISION 


143 

144 

200 



c 

5 ^^ 


IT, or, S1FW0,S2FWD,A1FWD,A2FWD,S1BACK,S28ACK,A1BACK,A2BACK,DA18AK, 402 
lDA2BAK,OELPl,DeLP2,DELQL,OELQ2,DELRl,DELR2,OELSl,DELS2,DELTl,DELT2 403 
2,DELC1,DELC2,CI,C2 

DOUBLE PRECISION YI,Y2,S1,S2,A1,A2,Y1P,Y2P,S1P,S2P,AIP,A2P, 

IY1C,Y2C, SlC,S2C,AlC,A2C,YlQ,Y2Q,SlQ,S2Q,AlQ,A2Q,YlT,y2T,SlT,S2T, 

2AIT,A2T, DYI,OY2,OS1,DS2,DA1,DA2,DYIP,DY2P,DS1P,DS2P,DA1P,DA2P, 
3DYiC,DY2C,DSlC,DS2C,DAlC,DA2C,DYlQ,DY2Q,DSlQ,DS2Q,DALQ,DA2Q, 

4DY1T, DY2T,DS1T,0S2T,DA1T,DA2T 

DOUBLE PRECISION EEfVV 

COMMON T ,DT 

COMMON SIFWD,S2FW0,AIFW0,A2FWD,S1BACK,S2BACK,A1BACK,A28ACK,DA1BAK, 402 
1DA2BAK,DELP1,DELP2,0ELQ1,0ELQ2,DELR1,DELR2,DELG1,DEL62,DELT1,DELT2 403 
2,DELC1,DELC2 

COMMON Cl,C2,AL,R,W,DDH,AA,AR,AHL ' 202 

EQUIVALENCE { Y1,T{ IH , {Y2,T(2) ), ( Si,T C3) ) , ( S2,T(4H , ( Al,T{5) ) ,{A2, 206 
IT(6)) ,(Y1P,T(7) V,(Y2P,T( 8) I,CSIP,T{9J ) , { S2P,T ( lOH , ( AlP.T C 11 H , ( A2 207 
2P,T(l2n,(YlC,T(13n,(Y2C,T(14)),(SlC,T(15n,{S2C,Ttl6H,CAlC,TC17 208 
3)), {A2C,T(18)) ,(YlQ,T(19n, {Y2Q,T(20) 1,(S1Q,T(21) ),(S2Q,T(22) ),( A1 
4Q,T(23)) , ( A2Q,T(24n, (YIT,T{ 25 ) I , I Y2T ,T{ 26) ) , ( SIT ,1(27) ) , (S2T ,T ( 28) ) , 

5) ) ,IA1T,T{29) ), IA2T,T(30) ) 

EQUIVALENCE ( DY1,DT { 1 ) ) , ( DY2,DT { 2 ) ) , { OSl , DT (3 ) ) , I DS2, DT< 4 ) ) , 

1{DA1,DT{ 5) ), (DA2,DT(6) ),{DYIP,DT(7)), (DY2P,DT( 8) ) , (DSIP,DT(9) ), 
2(DS2P,0T{10) ),(DA1P,DT{ ID), (DA2P,DT{ 12) ) , ( DY1C,DT( 13 ) ) , ( DY2C ,DT { 14 ) ) , 
34)),(0S1C,0T{ 15D,(DS2C,0TC16)),CDA1C,DT(17)),(DA2C,0T(18)), 

4{DY1Q,DT( 19) ),(DY2Q,DT(20)), (0S1Q,DT(21) ),tOS2Q.DT(22)) , 

5C0A1Q,0T(23) ),CDA2Q,DT{24) ), {DY1T,0T(25) ) , t DY2T,DT(26) ) , ( OSIT ,DTC27)), (OS 
67)),(DS2T,DT(28)),(0A1T,0T(29)),(DA2T,DT{30)) 


IF(I)52,52,54 

216 

FORWARD 54 

217 

VV(l,l)=Yl 

218 

VV{2,l)=Y2 

219 

VV(3, 1) = 0Y1 

220 

VV(4,1)=DY2 

221 

VV(5,1)=S1 

222 

VV(6,1)=S2 

223 

VV(7,1) = DS1 

224 

VV(8,1) = 0S2 

225 

VV(9,1)=A1 

226 

VV(10,1)=A2 

227 

VV(ll,l)=DAl 

228 

VV(12,1) =DA2 

229 

EEdt 1) = Y1P 

230 

EE(2, 1)=Y2P 

231 

EE(3,1)=DYIP 

232 

EE(4, 1)=DY2P 

233 

EE(5,1) = S1P 

234 

EE(6, iy=S2P 

235 

E£(7,1)=DS1P 

236 

EE(8,1)=DS2P 

237 

EE(9,1) = A1P 

238 

239 

EE{10,1)=A2P 

EEtll,l)=DAlP 

r 240 

241 

EE{12,l)=DA2P 

EE(1,3)= YIQ 

EE{2,3)= Y2Q 

EE(3,3)= DYIQ 



EEC4t3)= DY2Q 
EE(5,3H SIQ 
EE(6,3)= S2Q 
EE(7,3}= DSIQ 
eE{8,3)* DS2Q 
EEC9,3J= AIQ 
EE(10,3)= A2Q 
EEC 11,3)= OAIQ 
EEC 12,3)= DA2Q 
EEC1,5)=Y1C 
EE(2,5)=Y2C 
EEC3,5)=DY1C 
EEC4,5)=DY2C 
EEC5,5)=S1C 
EEC6,5)=S2C 
EE(7,5)=DS1C 
EE(8,5)=DS2C 
EE(9,5}=A1C 
EE(10,5)=A2C 
EE{li,5)=DAlC 
eE(12,5)=DA2C 
GO TO 56 
BACKWARD 52 
CONTINUE 

VVC 1,1) = Y1-VVC 1, 1) 

VVC2,1) = Y2-VV{2, 1) 

VVC3,1)=DY1-VVC3,1) 

VV{4,l)=0y2-VV{4,l) 

VVC5, 1) = S1-VVC5, 1) 

VV{6,1)=S2-VV(6,1) 

VVC7, 1)=DS1-VVC7,1) 

VV(8,1)=DS2-VVC8,1) 

VVC9, 1)=A1-VV(9,1) 

VVC10,1)=A2-VVC10,1) 

VVC11,1)=DA1-VVC11,1) 

VVC 12,1)=DA2-VVC 12,1) 

EeCl,5) = E6(l,5)-YlC 

EEC2,5)=EE(2,5)-Y2C 

EEC3,5)=EE<3,5)-DY1C 

EEC4,5)=EE(4,5)-DY2C 

EE(5,5)=EE(5,5)-S1C 

EE(6,5)=£E(6,5)-S2C 

EE17,5)=EEC7,5)-DS1C 

EE(8,5)=EEC8,5)-DS2C 

EE(9,5)=EEC9,5)-AIC 

EE(10,5)=EE(10,5)-A2C 

EEC ll,5)=EE( 11,5)-DA1C 

6E(12,5)=EE( 12,5)-DA2C 

EEC1,7)=-Y1P 

EEC2,7)=-Y2P 

EEC3,7)=-DY1P 

EEC4,7)=-DY2P 

EEC5,7)=-S1P 

EEC6,7)=-S2P 

EE(7,7)=-0S1P 

EEC8,7)=-DS2P 


259 

260 
261 
262 
263 

264 

265 

266 

267 

268 

269 

270 

271 
272 

273 

274 

275 
276 

277 

278 

279 

280 
281 

282 

283 

284 


2 



294 

295 

296 

297 

298 

299 

300 

301 

302 

303 

304 



EE(9,7)=-AIP 

EE(10,7)=-A2P 

EE(ll,7>=-DAlP 

EE112,7I=-DA2P 

EE(1,9>=-YIQ 

EE(2,9)=-Y2Q 

EE(3,9)=-DYiQ 

EE14,9)=-0Y2Q 

EE(5,9}=-SIQ 

EE{6,9)=-S2Q 

eE(7,9)=-0SlQ 

EE(8,9J=~DS2Q 

EE(9,9)=-A1Q 

EE(10,9)=-A2Q 

eE(ll,9J=-DAlQ 

EE(12,9)=-DA2Q 

EE{1,11)=-Y1T 

EE(2,11)=-Y2T 

EE(3,11>=-0Y1T 

EE{4,11}=-DY2T 

EE(5,11)=-SIT 

EE(6,1U=-S2T 

eE(7.1U=-0SlT 

EE(8, 11)=-DS2T 

EE(9,11}=-A1T 

EE(10,11}=-A2T 

EE( Ilfll)=-DA1T 

EE(12,lir=-DA2T 

EVEN COLUMNS 

DO 58 K=l»6 

DO 58 L=l,6 

EE(2*L-l,2*KJ=-EEt2*L,2*K-l> 

EE(2*L,2*K)=EE{2*L-lt2*K-ll 

CONTINUE 

Y1=0. 

DO 60 L=lt 12 
Yl=Yl+VViLf lJ*VV(Lf II 
CONTINUE 

IFCIPRINT.GT.IIGO TO 85 
PRINT B0,(VV(L,1I,L=1,12I 
F0RMATI/2X,6E16-8 /2X,6E16.8) 
CALL GAUSSIEEf 12,12tVV,l,l} 
Y2=0. 

DO 62 L=1,L2 
Y2=Y2+VVlLt 1V*VVCL,1) 

CONTINUE 

IFdPRINT.GT.DGO TO 90 
PRINT 80.{VV{Ltl)»L=l»12l 
CONTINUE 
0ELP1=VV( Itl) 

DELP2=VV(2,1} 

DELQ1=VVI3,1) 

DELQ2=VV(4,1) 

DELCl=VVt5tl) 

DELC2=VV{6,1) 

0eLRl=VV(7,l) 


305A 

3058 

306 

307 


327 

328 
329 

330 

331 

332 

333 

334 

335 

336 

337 
338 

338 

339 

341 

342 

343 

344 

345 

346 

347 


349 

350 

351 

352 


353 

354 

355 
356 

357 

358 

359 



D6LR2=VV{8,1) 

0EL61=VV(9, 1) 

OEt62=VV(10,l) 

OELTl=VVCll,l) 

0ELT2=VV(12»1) 

GO TO 56 
56 RETURN 
END 

£IBFTC GASDP NOPRNT 

SUBROUTINE GAUSS { A,NN,N, B,MM, M 1 

DIMENSION A(12,12),B(l2,l),IPIV0T{12r,PIV0T(l2»,IND£XC12#3J 
DOUBLE PRECISION A, 8 

EQUIVALENCE (IRQWtJRQWJ, { ICOiUM, JCOLUM ) , (AMAX, T, SWAP) 

C INITIALIZATION 

15 DO 20 J=ltN 
20 IPIV0T(J)=O 
30 DO 550 1=1, N 

C SEARCH FOR PIVOT ELEMENT 

40 AMAX=0.0 
45 DO 105 J=1,N 

50 IF ( IPIVOTC J)-l) 60, 105, 60 
60 DO 100 K=1,N 

70 IF( IPIV0TlK)-l)8O, 100,740 
80 IF{ABS{AMAX)-DABS(A(J,K)) 185,100,100 
85 IROW=J 
90 TCOLUM=K 
95 AMAX=A(J,K) 

100 CONTINUE 
105 CONTINUE 

110 IPIVOT(ICOLUM)=IPIVOT{ ICOLUMHl 
C INTERCHANGE ROWS TO PUT PIVOT ELEMENT ON DIAGONAL 

130 IFIIROW-ICOLUM) 140,260,140 
140 CONTINUE 
150 DO 200 L=1,N 
160 SWAP=A{ IR0W,L1 
170 A(IROW,L)=A(ICOLUM,L) 

200 A{ICOLUM,L)=SWAP 
205 IF(M) 260, 260, 210 
210 DO 250 L=l, M 
220 SWAP=B(IROW,L) 

230 81 IROW,L)=B{ICOLUM,L) 

250 B( ICOLUM,L INSWAP 

260 INDEXd, 1) = IR0H 

270 INDEK{I,21=IC0LUM 

310 PIV0T(I1=A( ICOLUM, ICOLUM) 

C DIVIDE PIVOT ROW BY PIVOT ELEMENT 

330 A( ICOLUM, IC0LUM)=1.0 
340 DO 350 L=1,N 

350 A( ICOLUM, LI =A( ICOLUM, L)/PIVOTCn 
355 IF(M) 380, 380, 360 
360 DO 370 L=1,M 

370 BI ICOLUM, L)=B{ ICOLUM, L)/PIVOT( I) 

C REDUCE NON-PIVOT ROWS 

380 DO 550 L1=1,N 
390 IFILl-ICOLUM 1400, 550, 400 
400 T=A{L1, ICOLUM) 


360 

361 

362 

363 

364 

365 

366 

367 


F4020007 

! 

F4020012 ! 

F4020013 

F4020014 

F40 20016 

F4020018 

F40200i9 

F4020020 

F4020021 ' 

F4020022 

F4020024 i 

F4020025 * 

F4020026 

F4020027 

F4020028 

F4020029 

F4020031 


F4020035 

F4020036 

F4020037 

F4020038 

F4020039 

F4020040 

F4020041 

F4020042 

F4020043 

F4020044 

F4020045 

F4020046 

F4020049 

F4020051 

F4020052 

F4020053 

F4020054 

F4020055 

F4020056 

F4020058 

F4020060 

F4020061 

F4020062 

, ' ' I 



420 A(L1, ICOLUM)=0.0 
430 00 ^50 L=ltN 

450 A(L1,L)=A(L1,L)-A{ IC0LUM,LJ»T 
455 IFCM) 550, 550, 460 
460 DO 500 L=1,H 

500 B(Ll,L)=B(Ll,L)-BtICOLUM,L}»T 
550 CONTINUE 
C INTERCHANGE COLUMNS 
600 DO 710 1=1, N 
610 L=N+1-I 

620 IF lINDEX(L,l)>INDEXtL,2) J 630, 710, 630 
630 JR0W=IN0EX(L,1} 

640 JC0LUH=INDEX(L,2I 
650 DO 705 K=1,N 
660 SWAP=ACK, JROWl 
670 A{K,JRO«)=A(K, JCOLUM) 

700 A (K, JCOLUM )=SWAP 
705 CONTINUE 
710 CONTINUE 
740 RETURN 
END 

£ENTRY 

I 24 30 1 15 15 1 

-0025 .05 1. .5 

.1000 .1000 1.0000 2 

0.0 l.O O.lOQO 

1.3000000015.00000000 4500.0 

-0.3098740794955324D 01 -0. 2234099599445573D-01 
-0.5765141780835510D 01 0.257235041956644 ID 02 

-0.4662111630903466D 04 0.15995141637509980 05 

-0.6787564833112370D 06 0.35701167022128940 00 


F4020063 
F4020064 
F4020065 
F4020066 
F4020067 
F4020068 
F4020069 
F4020071 ! 
F4020073 1 
F40 20074 i 
F4020075 1 
F4020076 * 

F4020078 

F4020079 

F4020080 

F4020081 

F4020082 

F4020083 

F4020084 


-0.78106392564845530 01 
-0.1058417864732388D 02 
-0.99567153004784980 04 
0.6071817 102296290D-02 



